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FORWARD 


This  is  the  first  part  of  a two  part  report  dealing  with  Techniques 
and  Applications  of  Computer-Aided  Circuit  Simulation.  The  topics 
discussed  in  Part  I are  the  backbone  of  a seminar  course,  EE  417-2,  which 
was  given  Autumn  1973.  Support  for  the  seminar  from  the  General  Motors 
Research  Laboratories  is  gratefully  acknowledged. 

During  the  course  of  EE  417-2  the  fundamentals 
developed.  At  the  same  time  eight  doctoral  students  used  CAD  program 
SPICE  to  investigate  specific  modeling  and  circuit  design  problems  which 
are  important  in  their  thesis  studies.  These  specific  problems  are 
coupled  with  projects  supported  by  the  NIH,  JSEP  and  ONR. 

A primary  result  of  the  seminar  has  been  to  bring  together  theoretical 
and  experimental  information  on  CAD  and  its  application  which  has  thus  far 
been  difficult  to  obtain. 

Part  I is  a concise  statement  of  CAD  techniques  with  examples  pertin- 
ent to  SPICE.  Part  II  reports  the  results  of  the  investigations  by  the 
doctoral  students  Involved  in  the  seminar.  The  titles  of  the  individual 
projects  are  given  in  the  last  subsection  of  Part  I. 

The  seminar  participants  are  gratefully  acknowledged  for  their 
individual  contributions  to  this  report  as  well  as  for  the  collective 
feedback  provided  by  the  seminar.  Special  thanks  goes  to  Mr.  Peter  Slapnicar 
and  Mr.  Tak  Young  who  were  major  contributors  to  the  seminar  via  discussions 
of  CAD  techniques  and  program  developments. 
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INTRODUCTION 

The  need  for  computer  simulation  is  steadily  increasing  in  the  design  of 
Integrated  Circuits  and  systems  using  IC’s.  The  process  dependence  of  component 
parameters  and  parasitic  device  effects  cannot  be  adequately  breadboarded  to 
test  performance  for  IC's.  Computer-Aided  Circuit  Analysis  provides  an  efficient 
and  flexible  means  to  evaluate  the  performance  of  circuit  designs  prior  to  the 
costly  and  time-consuming  fabrication  process.  Component  values  and  process- 
dependent  device  parameters  can  be  varied  to  optimize  circuit  performance  and 
improve  fabrication  yields.  However,  to  effectively  use  these  computer  tools, 
the  user  must  be  able  to  supply  the  programs  with  realistic  models  and  model 
parameters.  Furthermore,  a basic  understanding  of  the  techniques  used  for 
computer  circuit  simulation  can  allow  the  designer  to  critically  assess  the 
validity  of  the  simulated  results. 

• 

This  report  presents  material  and  examples  pertinent  to  the  enlightened  use 
of  sophisticated  computer  circuit  analysis  programs  for  the  design  of  Integrated 
electronic  circuits  and  systems.  The  examples  chosen  for  this  discussion  are 
developed  with  program  SPICE.  SPICE  is  a nodal  analysis  program  which  offers 
nonlinear  dc,  nonlinear  transient  and  linear  ac  analysis  of  electronic  circuits 
in  a single  program.  Free  format  input,  built-in  device  models  (diodes  and  bipolar 
and  field-effect  transistors)  and  circuit  nesting  features  make  the  program  easy 
to  use.  The  sparse  matrix  structures  and  the  use  of  implicit  integration  and 
adjoint  network  techniques  help  to  make  SPICE  an  efficient  simulation  tool. 

Fundamentals  of  the  SPICE  program  structure  will  be  discussed  in  the  sub- 
sections which  follow.  The  basic  analysis  approach,  data  structure  considerations 
and  a method  for  numerical  interaction  will  be  presented.  Next,  the  method  of 
adjoint  network  calculations  will  be  developed  and  applications  will  be  described. 

Finally,  and  most  importantly,  advanced  MOS  and  bipolar  transistor  models  will  be 
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be  developed  as  they  are  implemented  in  SPICE.  The  need  for  accurate  models  is 
of  paramount  importance  and  the  requirements  to  correctly  determine  parameters 
for  theae  models  will  become  clear  as  one  reads  the  project  reports  in  subsequent 
sections . 
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COMPUTER-AIDED  CIRCUIT  ANALYSIS  OVERVIEW 
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To  initiate  the  discussion,  the  circuit  example  shown  in  Figure  1 will  be 
used  as  a vehicle  to  illustrate  the  techniques  of  computer  aided  circuit  analysis. 

The  three  input  waveforms  shown  in  Figure  1 drive  the  circuit  so  that: 

a)  the  dc  steady-state  bias-point  is  achieved, 

b)  the  small  signal  steady  state  behavior  is  determined  (as  a function  of 
f requency ) , and 

c)  the  transient  output  waveform  is  calculated. 

Each  of  these  analysis  modes  requires  the  solution  of  the  nonlinear  set  of  equations 
associated  with  the  transistor.  For  cases  a)  and  b)  the  solution  with  nonlinear 
elements  must  be  found  only  once.  For  the  small  signal  analysis  (case  b))  the 
subsequent  frequency  calculations  use  the  model  values  linearised  about  the  bias 
point.  However,  for  case  c),  the  transient  analysis,  the  set  of  nonlinear  equations 
* must  be  solved  at  each  time— point  in  the  analysis.  The  solution  of  the  nonlinear 

equations  associated  with  the  device  model  is  thus  the  backbone  of  computer  circuit 

K- 

-.v'j  simulation.  Efficiency  in  achieving  this  solution  is  a major  concern. 

The  nonlinear  dc  solution  Is  achieved  by  iterative  solution  of  the  linearized 

equivalent  circuits.  The  so-called  "linearized  equivalent  circuit"  is  constructed 

using  the  first-order  terms  of  the  Taylor  series  expansion  of  the  nonlinearities 

about  some  initial  point.  This  particular  iterative  procedure  - often  called  Newton- 

(2) 

Raphson  iteration  - solves  for  the  unknown  voltages  in  the  Taylor-series  expansion, 
updates  the  expansion  about  the  newly  found  value  and  continues  the  process.* 

* 

Reference  (2)  is  included  in  Appendix  I to  provide  further  comments  and  references 
on  techniques  of  Computer-Aided  Circuit  Analysis. 


_____ 





Each  linearization  is  made  using  the  solutions  from  the  previous  iteration.  The 
iterations  continue  until  the  solution  converges  to  within  some  specified  tolerance 
of  the  value  from  the  previous  iteration  - for  example  lO^V  at  every  node.  All 
energy  storage  elements  are  replaced  by  their  dc  equivalents  for  the  dc 

solution.  The  network  equations  themselves  (with  linearized  matrix  entries  for 

(2) 

the  model)  are  solved  using  a modified  form  of  Gaussian  elimination.  For  SPICE, 
the  equations  are  written  based  on  Kirchhoff's  Current  Law  for  all  the  nodes  in 
the  circuit  (l.e.  nodal  analysis).  The  simplest  Gaussian  elimination  consists  of 
constructing  an  upper  triangular  matrix  from  the  original  nodal  admittance  matrix, 
Y,  and  determining  the  node  voltages,  in  sequence,  starting  from  the  bottom  of  the 
triangular  matrix.  This  final  solution  method  is  called  backward  substitution. 

The  SPICE  implementation  uses  a matrix  transformation  to  achieve  two  matrices  u and 
L which  are  upper  and  low  triangular  with  the  property  that  LU  ■ Y.  The  voltage 
vector  is  solved  for  by  simple  matrix  operations  involving  L 1 and  U 1 . The  matrix 
operations  are  easy  owing  to  the  triangular  structure  of  the  matrices . 

The  procedures  described  above  are  the  basis  of  essentially  all  nonlinear 
computer-aided  circuit  analysis.  For  small  signal  analysis,  with  frequency  as  a 
variable,  the  need  for  Iterations  is  eliminated  once  the  nonlinear  dc  solution  is 
found.  Specifically,  the  model  is  linearized  about  the  dc  operating  point.  The 
ac  analysis  involves  arithmetic  operations  with  complex  arguments.  However,  the 
basic  approach  of  the  L-V  tranai ormation  is  used  to  solve  for  the  complex  voltage 
vector.  The  energy  storage  elements  have  been  added  to  the  Y matrix  as  complex 
admittances • 

For  transient  analysis  the  procedure  described  for  the  nonlinear  dc  bias 
solution  is  used  repetitively.  Since  energy  storage  elements  are  now  included  in 
the  circuit,  the  solution  of  the  integral-differential  equations  must  be  found. 
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The  solution  of  the  integral-differential  equations  proceed  using  numerical 

(2) 

integration  techniques.  For  example,  the  branch  relationship  for  an  inductor, 

in  integral  form,  is: 

t+A 

\ (t+A)  * \ (t)  + if  vl(t)  dT 

and  the  Integral  can  be  approximated  as: 

i^  (t+A)  - iL  (t)  + ^ VL  (t+A) 

The  choice  to  approximate  the  integral  using  the  yet  to  be  determined  voltage  at 

t+A  is  a form  of  implicit  integration.  Now  the  branch  relation  for  the  inductor 

at  time  t+A  has  a form  equivalent  to  a current  source,  i (t),  (a  constant  since 

L 

we  have  just  solved  for  its  value)  in  parallel  with  a conductance  of  value  L. 

A similar  branch  relation  exists  for  capacitive  elements.  The  result  is  that  at 
each  new  time  point  t+A,  the  energy  storage  elements  are  replaced  with  their  equiva- 
lent circuits  which  approximate  the  integral-differential  branch  relationships. 

The  nonlinear  solution  of  the  circuit  equations  can  then  proceed  at  that  time  point, 
Just  as  originally  described-  for  the  nonlinear  dc  condition. 

Figure  1 and  the  above  discussion  have  been  used  only  to  introduce  the  elements 
of  computer-aided  circuit  analysis.  The  following  subsections  will  expand  these 
ideas  and  illustrate  their  implementation  using  the  example  in  Figure  1 and  the 
algorithms  in  program  SPICE. 
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MODEL  LINEARIZATION 
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Modeal  linearization  is  the  first  step  in  obtaining  the  nonlinear  dc  bias 
point  for  an  electronic  circuit  (for  example  the  circuit  shown  in  Figure  1).  To 

show  more  precisely  how  the  linearization  proceeds,  let  us  begin  by  considering 

(3) 

the  bipolar  transistor  model  showin  In  Figure  2.  The  equations  associated  with 
the  elements  are: 


I _ ■ I («qVBE/  -1) 
CE  S kT 


*00  " Js  (eqVBC/w  -1) 


qICE 

C„  * G (V  ) + ■ T_ 

E JE  BE  kT  F 


CC  * CJC(VBC)  + TF  • T* 


(1) 


(2) 


(3) 


(4) 


where  I is  the  base  transport  current  (I„  * Of  I * at  I„„  for  conventional  Ebers- 
S S F ES  R CS 

Moll  notation);  V and  V are  the  junction  voltages;  C and  C _ are  the  voltage 
BE  BC  JE  jC 

dependent  depletion  layer  capacitances  and  r and  r are  the  forward  and  reverse 

F R 

transit  times  which  describe  the  diffusion  capacitances  due  to  minority  charge 
storage . 

The  two  diodes  labeled  as  I and  I in  Figure  2 represent  the  forward 

CE/j3p  CC/p^ 

and  reverse  contributions  to  base  current,  using  the  exponential  relationships  given 

by  equations  (1)  and  (2).  The  I and  I terms  are  the  forward  and  reverse  base 

CE  CC 

transport  currents.  The  linearization  of  the  four  exponentially  nonlinear  elements 

shown  in  Figure  2 is  essential  to  initiate  the  iterative  solution  for  the  dc  bias 

point  of  the  circuit  shown  in  Figure  1.  To  illustrate  the  linearization  procedure, 

consider  the  exponential  diode  relationship  for  I as  is  shown  in  Figure  3a. 

CE/  a 

F 


w 


If  an  initial  operating  point  V is  chosen,  a linearization  about  that  point 

B£U 

can  be  made  by  the  parallel  combination  of  the  current  source  and  conductance 
shown  in  Figure  3b.  The  element  values  are: 


In  equation  (6)  the  first  term  represents  the  actual  base  current  for  a 

bias.  Performing  the  same  procedure  for  the  other  exponential  relationships, 

the  composite  linearized  dc  model  is  shown  in  Figure  4.  The  "0"  subscripts  have 

been  dropped  since  we  will  no  longer  restrict  ourselves  to  the  initial  values 

of  the  linearization.  For  the  linearization  of  I and  I it  should  be  noted 

05  CC 

that  the  conductances  are  in  fact  transconductance  generator  elements.  The 
importance  of  this  result  will  become  apparent  as- the  nodal  admittance  matrix 
is  generated  for  the  circuit  shown  in  Figure  1.  The  additional  terms  in  Figure 
4 have  values  given  by: 


G 


PIR 


' (7) 


rBCR 


G . V 
PIR  BC 


(8) 


(9) 
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The  two  current  generators  and  could  easily  have  been  expressed 

F BEF  R BCR 

in  a form  more  like  equation  (6)  by  defining: 

(11) 

(12) 

It  follows  immediately  that  * I __  and  3_ID_  » I , thus  only  the  form 

F BEF  CEF  R BCR  CCR 

given  in  Figure  4 is  used. 

The  means  for  applying  the  results  shown  in  Figure  4 can  now  be  summarized. 

Before  each  new  iteration  to  solve  for  the  dc  bias  point  of  the  circuit  shown  in 
Figure  1,  the  nonlinear  equations  describing  the  transistor  model  are  linearized 
about  the  voltages  from  the  previous  iteration  (or  some  default  values  for  the 
first  Iteration).  Figure  4 shows  the  linear  elements  used  to  represent  the  Bipolar 
Transistor.  In  the  next  subsection  the  entering  of  these  model  element  values  into 
the  y matrix  will  be  demonstrated. 

I 


ICEF*  ICE  ' °WF  ‘ VBE 
ICCR"  :CC  ~ GMR  • VBC 
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NODAL  CIRCUIT  ANALYSIS 


A.  MATRIX  FORMULATION 

Having  considered  the  linearization  of  the  bipolar  junction  transistor  model 
(as  an  example  of  a nonlinear  element  to  be  utilized  in  a circuit  design)  we  can 
now  proceed  to  see  how  the  linear  elements  are  entered  into  a formulation  which 
can  be  solved  efficiently  using  numerical  techniques.  The  approach  used  in  SPICE 
is  to  construct  and  solve  the  Kirchhoff  Current  Law  equations  for  each  node  of  the 
circuit.  The  circuit  shown  in  Figure  5 represents  the  bipolar  circuit  from  Figure  1 
with  the  transistor  replaced  with  its  linearized  equivalent  model  shown  in  Figure  4. 
All  energy  storage  elements  have  been  removed  from  the  circuit  for  the  dc  analysis. 
All  node  voltages  are  referenced  to  the  ground  node  (which  is  numbered  "0”).  The 
current  law  equations  for  nodes  1-5  are  written  below: 


i- 

— 

— — 

— — 

Node 

1: 

(“°IN'' 

0 

0 

0 

V1 

Iv 

IN 

Node 
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0 
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(-V 

(W 
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V4 

0 

Node 

5: 

0 

0 

0 

("V 

(C  ) 

L J 

la. 

1 

> 

H 

The  currents  are  taken  as  being  positive  out  of  each  node.  The  terms  I and  I 

V Ivi  V LL 

are  the  currents  through  the  voltage  sources.  These  current  values  are  determined 
by  other  circuit  elements  and  since  the  voltages  are  constant,  the  nodes  can  be 
rearranged  so  as  to  move  them  to  the  bottom  of  the  matrix.  The  five  equations 
were  written  directly  from  Figure  5.  The  following  rules  can  be  applied  to 
write  the  equations  knowing  only  the  nodal  connections  of  the  elements.  For  con- 
ductance between  nodes  i and  J the  value  is  added  to  the  diagonal  entries  y^  and 


and  is  subtracted  from  the  off-diagonal  entries  y^  and  y^.  For  a transconductance 
the  reference  polarities  and  node  numbering  shown  in  Figure  6 can  be  considered.  The 
transconductance  term  g^  is  added  to  the  matrix  terms  yik  and  y ^ and  is  substracted 
from  the  y^  and  y ^ terms.  This  usually  results  in  a non-symmetric  matrix  structure. 

For  the  case  shown  in  Figure  5,  i-J-0  for  the  G^  generator,  thus  only  the  y ^ term 
is  non-zero.  For  the  G^  generator  i-0  so  that  only  the  y^  and  y^2  entries  are  non- 
zero. It  follows  that  equation  (13)  could  have  been  written  directly  from  a list  of 
element  node  connections  and  prescribed  linearized  forms  for  the  nonlinear  device 
models.  However,  equation  (13)  is  not  necessarily  in  the  most  convenient  order  for 
efficient  solution  of  the  network  equations.  The  voltages  at  nodes  1 and  5 are  known, 
and  it  is  convenient  to  reorder  the  matrix  to  utilize  this  fact.  The  Y matrices 
encountered  in  nodal  circuit  analysis  generally  contain  many  zero  terms.  These  zero 
terms  can  easily  account  for  50-80Z  of  the  Y matrix.  It  is  possible  to  take  advantage 
of  the  sparsity  of  non-zero  terms  to  reduce  the  memory  required  to  store  the  Y matrix 

and  to  reduce  the  number  of  numerical  operations  required  to  solve  the  nodal  equations. 

(4) 

Appendix  II  contains  a detailed  description  of  one  particular  implementation  of 
what  has  come  to  be  called  "Sparse  Matrix  Techniques".  The  following  section  illustrates 
the  features  of  the  sparse  matrix  technique  for  efficient  storage  and  L-U  factorization. 

(4) 

The  circuit  from  Figure  5 is  used  as  an  example  and  the  notation  is  that  of  Berry 
B.  SPARSE  MATRIX  SOLUTION* 

The  first  step  is  to  reorder  the  rows  and  columns  (as  is  shown  in  equation  (14)) 
to:  1)  move  the  known  node  voltages  and  to  the  bottom  of  the  matrix;  and  2)  move 
to  position  shown  in  equation  (14^) . The  second  step  was  needed  to  minimize  fill-in 
during  the  L-U  transformation.  The  details  of  checking  for  fill-in  are  discussed  in 
the  Appendix.  The  resulting  equations  are: 


While  the  sparse  matrix  approach  can  offer  computational  advantages  it  can  also 
obscure  the  basic  solution  approach.  Because  of  this  conceptual  draw-back  it  is 
suggested  that  the  reader  skip  to  equation  (17)  for  the  first  pass. 
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The  3x3  matrix  shown  in  the  upper  left  corner  of  equation  (14)  is  the  only  portion 
of  interest  since  and  are  known.  What  is  done  next  is  to  create  pointers  to 
the  non-zero  entries  in  this  3x3  matrix.  This  is  best  illustrated  by  redrawing 
the  matrix  with  1 and  0 entries  corresponding  to  non-zero  and  zero  entries  respectively. 
Corresponding  to  columns  and  rows  1-3  of  equation  (14)  the  following  structure  exists: 


A column  array  is  created  which  contains  the  diagonal  Y matrix  entries  and  the  non- 
zero upper  and  lower  triangular  matrix  entries  in  a prescribed  sequence.  Two  arrays 
of  pointers  are  created  to  identify  the  matrix  indices  of  the  non-zero  terms  stored 
in  the  column  array.  This  is  done  in  two  passes:  First  pointers  to  every  non-zero 
entry  in  the  entire  matrix  are  created,  next  this  is  collapsed  to  a set  which  points 
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either  to  Upper  or  lower  triangular  entries.  In  the  "first  pass"  array  IUR 
has  n entries  corresponding  to  the  n-nodes  of  the  circuit.  The  ith  entry  in 
IUR  indicates  the  first  position  in  the  second  array,  IUC,  which  contains  a non- 
zero element  in  the  ith  row  of  the  Y matrix.  The  IUC  array  contains  the  column 
positions  of  the  non-zero  off-diagonal  entries  in  each  row  of  the  matrix  in  sequence, 
for  the  ith  row,  IUR  points  to  the  first  entry  in  IUC.  Further  entries  in  IUC 
contain  subsequent  column  numbers  of  matrix  positions.  For  the  3x3  matrix 
given  in  equation  (15),  the  IUR  and  IUC  arrays  contain  the  following  entries  on 
the  first  pass: 


IUR  (1)  * 1 


IUC  (1)  - 3 


IUR  (2)  * 2 


IUC  (2)  » 3 


IUR  (3)  « 3 


IUC  (3)  - 1 


IUC  (4)  . 2 


For  example,  IUR  (3)  "points"  to  IUC  (3).  The  IUC  (3)  and  IUC  (4)  entries  indicate 

pd 

that  columns  1 and  2 have  r.on-zero  terms  in  the  3 row.  The  "second  pass"  compares 
entries  corresponding  to  the  upper  and  lower  triangular  portions  of  the  matrix. 

A set  of  pointers  is  then  choosen  so  that  the  matrix  is  assumed  symmetric.  This 
new  set  of  pointers  correctly  selects  all  non-zero  entries  but  In  addition  may  point 
to  a limited  number  of  zero  terms  corresponding  to  the  non-syrametrical  entries.  'For 
the  example  given  by  equation  (15)  the  new  pointers  are: 


IUR  (1)  » 1 
IUR  (2)  « 2 


IUC  (1)  - 3 
IUC  (2)  a 3 
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Notice  that  only  n-1  rows  are  now  needed.  Considering  the  upper  portion  of  the  matrix, 
IUR  (2)  points  to  IUC  (2)  which  in  turn  indicates  that  the  3rd  column  of  row  two 
in  upper  triangle  is  non-zero.  The  same  pointers  refer  to  the  lower  triangle  in 
that  the  IUR  entries  now  go  column  by  column  and  the  IUC  entries  tell  the  row  of 
non-zero  entries.  Thus  IUR  (2)  indicates  that  in  the  2nd  column  we  must  go  to 
IUC  (2)  to  find  non-zero  entries  in  the  lower  triangle  of  Y.  IUC  (2)  tells  us  that 
the  entry  in  the  3rd  row  of  column  two  is  non-zero.  To  see  how  the  pointers  and 
column  array  fit  together,  let  us  enumerate  the  Y matrix  entries  in  equation  (13). 


The  overall  data  structure  is: 


IUR  (1)  ■ 
IUR  (2)  » 


pointers 


col urn  arrays 
*11  * A(1)  ) 


y m A(2)  ) diagonal  entries 

22  ( 


*IUC  (1)  - 3 *yi3  . A(4)  1 

*IUC  (2)  - 3v-.w *y23  - A(5)  ) 


y33  -A(3> 

*yi3  " A(4) 


>y23  - A(5)  • 
-A(6)  | 

N^32  -A(7)  > 


upper  triangle  entries 


lower  triangle  entries 


The  diagonal  entries  are  of  known  number  and  location.  The  index  of  IUC  tells 
the  total  number  of  entries  in  either  triangle.  Thus  the  starting  location  of  the 
upper  triangle  in  the  column  array  is  n+1.  For  the  first  lower  triangle  entry,  the 
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position  is  at  n+1  incremented  by  the  maximum  index  for  IUC. 

This  rather  long  procedure  has  given  us  a sparse  pointer  structure  which  locates 
all  non-zero  entries  in  the  Y matrix.  Using  these  pointers,  the  following  algorithm 
can  be  used  to  calculate  the  L and  U matrices  such  that  LTJ  ■ Y.  The  procedure  must 
be  executed  n times  (where  n is  the  number  of  circuit  nodes)  and  the  superscript 
indicates  the  ith  or  i-lst  row-column  step  in  the  procedure. 


U1J  " yi(j_1)  f°r  811 


y . (i-i) 

1 , « J1  for  j>l 

J1  u 

ii 


and  y^-  y^"1’  - lJt  ulk  for  J and  lc  > i 


Once  the  entries  in  the  L and  U matrices  are  known,  the  solution  for  node 


voltages  is  straightforward.  Kamel y, 


Yv  - LIV  - i 


Inverting  L and  multiplying  both  sides  of  (17)  by  L , one  obtains: 


Vv  - l"1!  4 i*-  (18) 

a 

Inverting  U and  multiplying  both  sides  of  (18)  by  U_1,  the  desired  voltage  vector 
is  obtained: 


v ■ U_1i*  ■ U-1L-1i  (19) 

Thia  procedure  is  used  for  nonlinear  dc  and  transient  analysis  with  Iterative 
solutions  as  models  are  updated  at  each  new  set  of  node  voltages.  The  procedure 
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is  also  appropriate  for  ac  analysis  using  complex  arithmetic  and  linearized  model 
entries. 


The  circuit  example  given  in  Figure  5 is  rather  simple  and  does  not  give  a fair 
picture  of  the  storage  and  computation  savings  that  are  possible  using  the  complete 
Berry  algorithm.  A more  typical  circuit  to  be  analyzed  and  the  resulting  matrix 
structure  are  shown  in  Figures  7a)  and  b) . For  this  case  the  number  of  zero-valued 
entries  is  substantial.  The  X's  in  Figure  7b)  indicate  the  original  Y matrix  entries 
corresponding  to  the  uncircled  node  numbers.  The  0's  indicate  "fill-in"  in  che  L and 
U matrices  resulting  from  the  decomposition.  Figure  7c)  shows  the  matrix  structure 
for  the  same  circuit  but  with  the  nodes  renumbered  as  shown  by  the  circled  node  numbers. 
It  is  clear  from  Figure  7c)  that  fill-in  is  considerably  reduced  as  is  shown  by  the 
substantial  decrease  in  "0"  entries.  For  the  example  cited  in  Figure  7,  the  following 
reductions  in  arithmetic  operation  counts  can  be  achieved: 


DIV 

MULT  and  ADDS 

Straight  LU  decomposition 

378 

6930  ea 

Renumbering  and  no  sparse  pointers 

134 

746  ea 

Renumbering  and  sparse  pointers 

63 

169  ea 

The  reordering  aspect  of  che  Berry  algorithms  is  described  in  Appendix  II. 

' The  primary  purpose  of  the  discussion  of  sparse  matrix  techniques  has  been  to 

illustrate  the  creation  of  the  pointer  structure  needed  for  sparse  storage  and  the 
steps  needed  for  che  L-U  factorization.  A general  result  which  can  be  stated 
regarding  the  computational  efficiency  of  using  Sparse  Matrix  Structures  is  that  che 

3 

long  operations  count  (multiplications  for  example)  descreses  from  n/^  to  a number 
more  nearly  proportional  to  n,  where  n is  the  number  of  nodes  in  the  circuit. 
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NUMERICAL  INTEGRATION 

As  was  suggested  in  the  overview,  numerical  integration  approximations  are 
used  to  solve  the  integral-differential  circuit  equations  in  the  time  domain. 

The  form  of  the  integration  rule  which  is  used  can  make  a significant  difference 

in  both  the  speed  and  accuracy  of  the  computer-simulated  results.  If  we  consider 

(2) 

the  integration  of  a function  f(t),  three  common  approximations  which  can  be 
made  are: 

A . f(t) 

A . f (t+A)  (20) 

| [f(t>  + f(t+A)j 

These  integration  rules  are  celled  Forwsrd-Euler  (F-£),  Backward-Euler  (B-E)  and 

Trapezoidal  (TR)  respectively.  The  first  method  is  an  explicit  scheme  in  that  it 

uses  a known  function  value  to  approximate  the  integral.  The  Backward-Euler  and 

Trapezoidal  are  implicit  in  that  they  use  a not-yet-determined  value  of  the  function 

to  approximate  the  Integral.  The  use  of  the  implicit  integration  rules  has  been 

an  Important  advance  for  computer  circuit  simulation  owing  to  the  improved  stability 

(2) 

of  these  methods  over  Forward-Euler  for  a given  step  size.  It  is  well  known 
that  the  explicit  integration  may  proceed  no  faster  than  the  smallest  time  constant 
of  the  circuit  to  maintain  a stable  solution.  In  the  discussion  presented  here, 
the  accuracy  and  stability  properties  of  B-E  and  TR  integration  formulae  will  be 
developed  and  compared  with  that  for  the  F-E  method.  j 

First,  let  ua  understand  exactly  how  the  integration  formulae  are  applied 
computer-aided  circuit  analysis.  It  has  previously  been  shown  that  for  an  inductor 
the  branch  relationship  requires  that: 


t+A 


j f (r)dr 
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r t+A 


iL(t+A)  - it(t)  + i 


vL(r)dT 


Hence  our  numerical  approximations  to  the  Integral  can  be  applied  directly. 
The  results  are: 


+! 

(F-E) 

i (t+A)  « 

1 1 

( i (t)  + 7 v (t+A) 
j L la  -U 

(B-E) 

(21) 

1 

[ 1L(t>  + k VL(t)  + k VLCt+ ^ 

(TR> 

These  three  results  can  be  represented  as  the  branch  equivalent-circuits  shown 
in  Figure  8a.  These  element-values  are  entered  in  the  Y matrix  at  each  time 
point  and  the  solutions  for  t+A  is  found  using  Newton  iterations.  Note  that 
the  current  generator  values  may  change  as  time  proceeds.  If  the  program  includes 
variable  time- increments , A,  then  the  conductances  also  change  value  with  time. 

Constant-valued  capacitive  elements  impose  the  following  branch  relationship: 


V*>  - 4C<,) 


(22) 


which  can  be  rewritten  as: 


t+A 

u 


ic(r)dT  -vc(t+A)  - vc(t) 


(33) 


Again  we  can  apply  the  integration  rules  to  obtain  the  relationships: 
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1 


/ 

A 

vc(t) 

+ £ *-<t> 
c c 

(F-E) 

v (t+A)  * Jv  (t) 

+ A i (t+A) 

(B-E) 

c \ c 

c c 

(vc(t) 

+ A i (t ) + 

A i (t+A) 

(TR) 

2C  C 

2C  C 

These  results  can  be  represented  as  the  branch  equivalent  circuits  shown  in 
Figure  8b.  The  circuits  shown  in  parentheses  are  the  Norton  equivalents  which 
are  more  appropriate  for  nodal  analysis. 

Thus  far,  only  constant  valued  energy  storage  elements  have  been  considered. 
Assume  that  the  branch  relationship  is  now  defined  by  the  equation: 

ic(t)  M!tQ(vc(t))  (25) 

h sides  of  the  equation  yields 

t+A 

ic(r)dr  - Q(vc(t+A))'  - Q(vc(t))  ^ (26) 

t 


For  this  case  the  integral  must  be  approximated,  as  before,  but  in  addition  the 
nonlinear  Q relationships  must  be  linearized.  The  result  is  that  Newton-Raphson 
iterative  procedure  must  now  be  applied  to  nonlinear  energy  storage  elements  as 
well  as  to  the  dc  transistor  model  nonlinearization. 

A primary  result  of  the  preceeding  discussion  is  that  equivalent  circuits 
are  created  to  approximate  the  integral-differential  equations  associated  with 
energy  storage  elements.  These  equivalent  circuits  are  used  with  the  standard 
iterative  matrix  solution  methods  developed  for  nonlinear  dc  analysis  to  solve  for 


circuit  behavior  in  the  time  domain.  The  question  of  relative  accuracy  and 
stability  of  the  several  integration  methods  (F-E,  B-E  and  TR)  is  now  discussed 
in  the  context  of  the  simple  circuit  example  shown  in  Figure  9.  The  voltage 


vR  is  the  variable  of  interest.  It  is  apparent  from  the  circuit  that 


so  that 


vR(t)  ® E - vc(t) 


v (t+A)  - E - v (t+A) 


We  can  now  apply  the  integration  formulae  given  by  equation  (24)  to  approximate 
v^Ct+A).  That  is,  using  equation  (24)  in  equation  (28): 


v,  (t+A)  m E -v  (t)  - 
R C 


- i (t) 

c cy  J 

(F-E) 

(B-E) 

(29) 

* lc(t)  +lc  ic(t+A) 

(TR) 

It  is  true  that  for  any  time  t (or  t+A): 


left)  - i,(t> 


V‘> 


Thus  using  equations  (27)  and  (30)  in  equation  (29)  one  obtains: 


Yt+A>  - %<t) 


— v (t) 

RC  RV 

A v (t+A) 
RC  R 


^ V°  * 25c  <™> 
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These  equations  can  be  rearranged  to  give  the  following  recursive  relationships: 


/V*>  ('  --r) 
V'*4'  * iV*’  l-±j\ 


'1  + 2t 


where  twRC.  For  a step  voltage  input,  the  initial  value  of  vR  is  vR(o)  ■ EQ. 

By  applying  the  formulae  given  by  equation  (32)  repetitively  for  n equal  A’s  the 


result  is: 


vR(nA) 


nr 


’V*k 


The  exact  solution  for  this  circuit  and  step  stimulation  is: 


v (t)  m E e 
R o 


(EXACT) 


so  that  for  t « nA: 


v (nA)  m E e T 
R o 


(EXACT) 


Several  comparisons  can  now  be  made  using  equations  (33)  and  (35).  For  a 
single  step  (nal)  the  accuracy  can  be  compared  for  several  values  of  u.  The  results 
for  A m ,1T  and  A * t are  given  below. 


ACCURACY: 

F.E. 

M) 

B.E. 

1 

1 

T 

— 

TR. 

EXA£T 

Y 

e 

A « O.lr 

0.  90000 

0. 90909 

0.  90476 

0.  90483 

A . T 

0.0 

0.5 

0.333 

0.368 

Clearly  it  can  be  concluded  that  the  trapezoidal  approximation  gives  the  most 
accurate  approximation  to  the  exact  solution  for  a given  time  step  A.  The  stability 
of  the  methods  can  be  demonstrated  in  an  approximate  sense  by  seeing  for  what  values 
A the  solution  begins  to  oscillate  with  odd  and  even  values  of  n.  The  results  given 
below  show  several  things . 


First,  the  solutions  exhibit  oscillations  only  for  F-E  and  TR  when  A is  greater  than 
r and  2t  respectively.  The  B-E  formulation  is  stable  in  the  sense  that  the  / l \ 


I + * 


I — 


recursive  expression  never  oscillates  in  sign.  The  second  comment  to  be  made  is 
with  regard  to  stability  £or  very  large  values  of  A.  This  consideration  shows 
that  F-E  not  only  oscillates  but  diverges  from  the  correct  answer.  Both  B-E  and 
TR  are  stable.  However,  for  TR  the  oscillations  still  persist  for  conditions 
such  that  A>2t. 

Program  SPICE  uses  the  trapezoidal  integration  rule.  Thus  the  solutions 
are  of  the  greatest  accuracy  for  the  correct  choice  of  time  increments.  However, 
for  increments  which  are  too  large  ringing  can  be  observed  in  the  waveforms.  This 
is  most  often  eliminated  by  decreasing  the  internal  program  integration  increments 
while  maintaining  the  same  output  plot  increments.  It  might  be  thought  that  the 
backward  euler  formulation  would  be  a better  approach  owing  to  its  unconditional 
stability.  However,  one  major  drawback  should  be  stated.  There  is  no  easy  way  of 
telling  when  accuracy  is  degraded  due  to  excessive  time-step  increments.  With 
trapezoidal,  the  oscillation  problems  occur  at  about  the  step  increment  where  accuracy 
is  also  degraded.  Thus,  by  using  a proper  increment  step  control  both  the  accuracy 
and  stability  can  be  maintained  for  the  trapezoidal  method. 
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ADJOINT  NETWORK  ANALYSIS 

In  the  preceeaing  subsections  numerical  and  network  analysis  techniques 

have  been  described  which  are  the  basis  of  computer-aided  circuit  analysis. 

In  this  section  the  principles  and  application  of  the  adjoint  network  will 

be  described.  Although  the  topic  may  at  first  sound  theoretical  and  unrelated 

( 

to  practical  circuit  analysis,  in  reality  the  technique  is  easy  to  apply  and 
has  had  a major  impact  on  integrated  circuit  design  using  the  computer. 

.Specifically,  adjoint  network  techniques  are  used  to  provide  efficient  noise 

(5)  ^ 

analysis,  dc  and  ac  sensitivities  of  circuit  outputs  to  parameter  variations 

and  finally  "design  optimization".  One  hesitates  to  use  this  last  term  owing 
to  the  many  misconceptions  as  to  its  meaning.  In  fact,  the  adjoint  network 
can  be  used  along  with  the  standard  network  techniques  to  perform  two  types 
of  analysis  which  together  might  be  called  design  optimization.  Sensitivity 
analysis  is  used  to  determine  the  gradient  of  a circuit  output  with  respect 
to  the  circuit  parameters.  A performance  index  can  be  defined  and  the  devia- 
tion from  this  index  can  be  minimized  using  the  gradient  information  and  some 

algorlth  to  seek  the  minimum.  For  example,  the  Fletcher-Powell  search  algorithm 
(7) 

is  commonly  used.  This  approach  is  called  design  optimization.  Stated  more 
exactly,  one  can  say  that  sensitivity  analysis  is  used  iteratively  with  some 
search  algorithm  to  minimize  an  error  function.  A minimum  error  should  thus 
guarantee  an  "optimum"  design.  Unfortunately,  circuit  parameter  tolerances 
about  this  optimum  design  can  cause  a poor  production  yield.  This  is 
especially  true  for  integrated  circuits  where  exact  component  values  are 
difficult  to  guarantee.  For  IC's  it  may  be  more  advantageous  to  define  an 
optimum  design  in  terms  of  a production  yield.  Adjoint  network  analysis 
techniques  can  be  used  to  perform  worst-case  and  statistical  analysis.  In 
both  cases  the  adjoint  network  provides  an  efficient  means  to  analyze  the 
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effects  of  numerous  possible  circuit  variations  The  added  information 


& 


from  such  simulations  can  lead  to  an  optimum  component  value  and  tolerance 

assignment.  In  this  cas-e  the  objective  might  be  to  optmize  yield.  One 

can  see  that  the  efficient  analysis  of  many  possible  component  values  is 

* 

the  key  to  optimizing  a design.  Use  of  the  adjoint  network  provides  these 
needed  analyses  with  a minimum  increase  in  computational  expense. 

To  compute  the  sensitivities  for  a variable,  say  the  changes  in  voltage 
AVj  across  a current  source  1^,  we  assume  that  an  expression  for  can  be 
written: 


f(pl’P2  • • • pp> 


(36) 


where  the  p's  represent  the  parameters  which  give  rise  to  the  observed  AV^. 
The  change,  AV^  is  then  given  by: 


A V 


. ■ E t.  • 


(37) 


kal 


df 


and  each  term  -r—  represents  the  sensitivity  of  V to  the  change  in  parameter 
3pk  1 

p^.  Our  task  then  is  to  find  an  expression  of  the  form  of  equation  (37)  and 
to  identify  the  proper  sensitivity  terms.  Toward  this  end  the  adjoint  net- 
work is  now  introduced.  This  mathematical  formulation  along  with  the  proper 
reordering  and. Identification  of  terms  will  facilitate  the  desired  result. 

The  starting  point,  for  considering  the  adjoint  network  is  Tellegen's 
(8) 

Theorem.  The  theorem  states  that  if  we  have  two  topologically  identical 
network  (i.e.  the  same  graph  structure,  branch  numbers,  node  numbers  and 

A 

orientations)  - call  them  T|  and  1)  with  V^-I^  and  voltage-current 


It  should  also  be  realized  ,that  the  real  key  is  to  have  a good  design  to 
begin  with.  However  even  initial  designs  can  be  marked  improved  by  knowing 
sensitivity  information. 
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r 

B - 


identify  -tions  respectively  - then 


V,  * 0 
k k - 


k*l 


(38) 


b 


la  Vk  5 0 

kwl 

where  b is  the  total  number  of  branches  of  the  network.  Figure  10  shows 
the  networks  *1  and  ^ schematically,  along  with  a network  n-fAr|  which  represents 
incremental  pertubations  to  the  voltage  and  current  vectors  of  network  q. 

All  three  networks  satisfy  Tellegen's  theorem  and  it  can  be  easily  shown  by 
application  of  equation  (38)  that  Aq  also  satisfies  the  theorem  so  that: 


AV,*,  m 0 
k k - 


kml 

b 


(39) 


kml 

Subtracting  the  above  equations  from  one  another  it  follows  that: 


V Av,_«,_  - Y..AI 


La 

kml 


k k 


k k 


(40) 


This  equation  is  fundamental  to  determination  of  sensitivities  and  the 
definition  of  the  adjoint  network. 

Consider  now  the  possible  branch  relationships  defined  by  Figure  11 
for  current  sources,  voltage  sources,  resistors  and  conductances.  Equation  (40) 
can  be  rewritten  in  terms  of  this  notation  with  summations  taken  over  the 
respective  number  of  branches  for  each  element  type. 
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(41) 


Z K*v  - V'v]  ♦ £K4:  - V'l. 

V 1 

®r*r  - V’J  * Z K4s  ' V“o  1 ; 0 

G 

The  exact  relationship  for  the  change  In  voltage  across  a resistor,  for  example, 
is: 

V%  ■ (V“r)  (V“«)  <42) 

If  the  second  order  term  in  dR^dl^  is  neglected  then  it  follows  that: 

4\  * '»  Wi  <43> 

using  equation  (43)  In  the  summation  over  resistive  branches  In  equation  (41) 
the  result  Is: 

£ ("»«»  - V‘r)  ■ E K’r4r  * Vr4r  - V'rI  (44> 

H - R 

* S [*v«*»  ♦ My»  ~ M I 

(a)  (b)  (c) 

The  three  terms  identified  as  (a),  (b)  and  (c)  are  now  considered.  We  want 
to  determine  variations  based  on  (a),  however  term  (b)  also  enters  Into  the 
summation  and  is  an  unknown  quantity.  By  choosing  the  relationship  in  (c) 
so  that: 

’»  V.  <45> 

then  the  second  term  In  equation  (44)  goes  to  zero  and  only  the  desired  term 
Is  left.  What  has  been  done  la  to  define  the  adjoint  branch  relationship 
given  by  equation  (45)  so  that  this  happens.  Notice  that  this  so-called  adjoint 
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network  q satisfies  Tellegen's  theorem  with  regard  to  q but  the  branch 
relationships  may  be  different.  The  same  result  is  obtained  for  conductances 
so  that: 


LKW'o)  * £-“oVo 


and 


* • G-T 
G G G 


(46) 


(47) 


It  can  be  seen  from  equations  (45)  and  (47)  that  resistances  and  conductances 
translate  into  the  same  relationships  in  the  adjoint  network.  Using  equations 
(44)  and  (46)  in  equation  (41): 


Elv»-vvl 

V I 


* S va  * E -‘sv-  ■ 0 


G G 


which  can  be  rearranged  so  that: 


2 V:v  + 2 “Vi 

V I 


v&v + + Y -Ao-V-T-  (48) 

L V V Ls  II  ZrfR  R.R  G G G 

V I R G 

It  should  be  noted  that  all  summations  on  the  right-hand  side  of  equation  (48) 

contain  variations  in  the  element  values  themselves  (i.e.  AV  , AI  , AR  and 

V I R 
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while  the  variations  on  the  left-hand  side  of  the  equation  pertain 

to  either  voltages  across  current  sources  or  currents  through  voltage 

soucres  (AV^  and  AIV  respectively). 

It  is  now  a fairly  straightforward  matter  to  show  how  equation  (48) 

can  be  reduced  to  the  desired  result  given  by  equation  (37).  We  have  defined 

branch  relationships  in  the  adjoint  network.  By  an  appropriate  choice  of 

stimulations  of  the  adjoint,  Y.  and  * , the  respective  <6  , Y , « and  Y_ 

V I V I R G 

can  be  solved  for.  If  we  choose  all  Y^  and  to  be  zero  except  for  the 
single  current  source  for  which  we  want  to  determine  output  voltage  sensitivity 
we  obtain: 


-"i  ♦. 
n n 


all  other 

\ and 

are  zero 


I>W  + Z>Yi 

v I 


(49) 


* e4Vr*'i.  * s-Vo’g 

R 


By  setting  <t>j  m -1  the  form  of  equation  (49)  looks  very  much  like  that  of 
n 

equation  (37) , In  fact,  the  sensitivity  of  V to  changes  in  R for  example 

X R, 

3f  n k 

is  then  — - « I ® . To  determine  these  sensitivities  one  appropriately 

®*‘r  Rk  Rk 
k 

stimulates  the  adjoint  network  and  solves  for  the  Y and  I vectors.  Since  the 
original  V and  I vectors  are  known,  the  sensitivities  can  be  constructed. 

Several  factors  must  still  be  discussed  before  the  adjoint  network  method 
is  in  a form  suitable  for  application  in  computer  aided  circuit  analysis. 
Branch  relationships  for  all  elements  (i.e.  storage  elements  and  controlled 
sources)  must  be  defined.  In  addition,  the  approach  must  allow  calculation 
of  sensitivities  across  elements  other  than  the  restricted  case  for  voltage 
and  current  sources  which  is  suggested  by  equation  (48),  Finally,  a means 
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for  the  implementation  of  the  method  must  be  described. 

The  definition  of  all  elements  in  the  adjoint  is  shown  symbolically 
in  Figure  12.  For  purposes  of  nodal  circuit  analysis  the  transconductance 
element  is  the  single  most  important  relationship  for  active  circuits  which 
was  not  considered  in  equation  (48).  For  a voltage-controlled  current  source 
with  input  and  1^  m 0 and  output  • gnAT^  with  V2  arbitrary,  the  adjoint 
relationships  are: 

*2  * ° 

- gm»2  (50) 

F , arbitrary 

These  relationships  are  shown  in  Figure  12. 

The  application  of  the  adjoint  method  to  calculate  variations  across 
elements  other  than  sources  is  a straightforward  matter.  If  the  current 
through  some  arbitrary  element  is  desired,  then  a zero-valued  voltage  source 
is  entered  in  series  Into  the  network  and  the  Al^  can  be  calculated  as 
described  above.  Similarly  voltages  across  elements  can  be  calculated  by 
adding  zero-valued  current  sources  in  parallel  with  the  elements  for  which 
AVj  is  desired.  The  net  result  is  that  the  original  circuit  is  modified 
appropriately  so  that  the  desired  output  ports  can  be  stimulated  in  the 
adjoint  network  to  determine  sensitivities. 

The  method  of  computer  implementation  for  the  adjoint  is  found  to  be 

A 

straightforward.  To  compute  F and  f,  the  adjoint  matrix  Y is  needed.  However 
since  T and  r\  are  interreciprocal: 

$«YT  (51) 

thus; 

YTF  > f (52) 
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and  using  the  L-U  factorization: 


tlVr  . f (53) 

It  should  be  clear  that  once  L and  U have  been  constructed  for  the  original 
network,  the  aubsequent  adjoint  analysis  to  solve  for  T and  $ is  considerably 
less  time  consuming.  In  fact  to  a good  approximation,  if  the  original  solution 
for  V takes  a unit  of  time,  the  subsequent  solution  for  Y take  .1  units. 

In  conclusion,  the  adjoint  network  is  an  efficient  means  to  compute 
circuit  sensitivities.  The  adjoint  network  results  from  an  application  of 
Tellegen's  theorem  and  the  proper  definition  of  branch  relationships  in  that 
network.  The  adjoint  admittance  matrix  is  obtained  simply  since: 


A 

Y 


T T T 
Y1  . 0 L 


(54) 


T T 

The  operations  to  obtain  U and  L are  trivial  once  U and  L are  known.  Thus 
the  adjoint  calculations  require  a minimal  amount  of  additional  analysis  time. 
The  sensitivity  calculations  are  the  basis  for  computer  circuit  optimization 
using  gradient  search  methods.  Sensitivities  can  also  be  used  to  perform  noise, 
statistical  and  worst-case  analysis. 
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BIPOLAR  AND  MOS  TRANSISTOR  MODELING 


(3 ) 

Tbe  linearization  of  a simple  bipolar  transistor  modelv  ' was  used 
earlier  to  demonstrate  model  implementation  for  computer-aided  circuit 

(9) 

analysis.  In  this  subsection  the  specific  features  of  the  Gummel-Poon 

(10) 

bipolar  transistor  model  and  the  Shichman-Hodges  MOS  transistor  model 

will  be  discussed.  However,  to  understand  the  need  for  such  advanc'ed  model 

formulations  consider  again  the  simple  bipolar  transistor  model  shown  in 

Figure  2.  The  collector  current  characteristics  for  this  model  are  shown 

in  Figure  13a.  The  additional  characteristics  in  the  figure  show  the 

assumed  voltage  dependences  for  I and  I and  the  resulting  current  gain 

B C 

dependence  on  1^,  Finally  the  transit  time  is  assumed  constant  with  current 
level.  The  junction  space-charge  capacitances  have  the  familiar  voltage 
depe.  ce: 


C(V) 


(55) 


where  C is  the  zero-bias  capacitance,  \ is  the  Junction  potential,  m is  a 
JO 

coefficient  ranging  typically  froir,  1/2  to  1/3  and  V is  the  Junction  voltage 

(positive  for  forward  bias).  When  the  experimental  behavior  of  integrated 

bipolar  transistors  is  observed,  the  characteristics  are  more  nearly  those 

depicted  in  Figure  13b,  Namely,  the  I vs  V behavior  shows  a non-zero 

C CE 

slope  in  the  normal-active  region  which  is  the  device  output  conductance 

resulting  from  base-width  modulation.  In  addition,  current  gain  is  not 

constant  with  collector  current  which  results  from  the  low-level  recombination 

and  high-level  injection  effects  shown  by  the  I and  I curves  as  a function 

C B 

of  V . Finally  the  high  current  effects  also  increase  the  transit  time 
BE 

resulting  in  a decreased  cut-off  frequency.  The  measured  junction  space- 
charge  capacitance  can  still  be  adequately  modeled  by  equation  (55).  To 
achieve  model  behavior  for  the  bipolar  transistor  which  correctly  predicts 
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the  characteristics  shown  in  Figure  13b,  second  order  device  effects  must  be 


£ 


accounted  for. 

A.  THE  BIPOLAR  TRANSISTOR 

The  Gununel-Poon  model  incorporates  the  second  order  effects  described 
above  within  the  basic  model  formulation  shown  in' Figure  2.  The  transport 


current  is  modified  in  the  Gummel-Poon  model  by  dividing  by  the  normalized 


based  charge  term,  Q . Q is  the  total  majority  charge  stored  in  the  base 
B B 


region.  Base  width  modulation  tends  to  decrease  Q and  the  collector  current 

B 


thus  Increases.  High  level  injection  adds  majority  charge  to  the  base,  thus 


increasing  Q and  decreasing  I . The  implementation  of  the  Gummel-Poon  model 
B C 


in  SPICE  defines  Q in  terms  of  the  following  parameters: 
B 


V V 

«!-**»*  ‘V21 

A B 


(56) 


I 

s 

■ 

,1 

rs 

exp 

\ — r1 

Jkr 

exp 

i— / _i 

(57) 


Ql  4* 


yV + ^ 


(58) 


where  I_,  I,  and  V are  defined  in  Figure  13b  for  the  normal-active  case.  The 
S k A 


parameters  I and  V are  similarly  defined  from  the  reverse-active  characteristics, 
kr  B 


The  I Intercept  is  the  same  for  forward  and  reverse  operation  by  the  nature  of 

s 


(9) 


the  transport  model  formulation.  The  Q term  represents  majority  carrier 

2 


base  charge  that  la  added  via  injection.  The  Q term  represents  space-charge 


widening  effects  on  Q , Equation  (56)  assumes  that  changes  in  base  charge  are 

B 


linear  with  Junction  voltage  thereby  implying  constant  junction  capacitances. 


The  resulting  expression  for  I -I  (see  Figure  2)  is: 


'a1*  • [exp  (l?)  ' exp 


(59) 
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Base  current  in  the  SPICE  Gummel-Poon  model  consists  of  four  terms. 


The  diode  terms  shown  in  Figure  2 represent  the  ideal  components.  That  is, 

these  base  current  terms  have  an  exp|iJjpj  dependence  and  give  curves  parallel 

to  the  I and  I curves  for  forward  and  reverse  operation  respectively. 

C E 

The  iderl  components  are  described  in  terms  of  3 and  3„„>  the  maximum 

FM  RM 

forward  and  reverse  current  gains.  The  non-ideal  components  represent  space- 

charge,  surface  and  other  recombination  effects.  Their  voltage  dependence  is  of  the 

form  exp  where  ne  is  the  parameter  for  the  base-emitter  junction  and 

nc  is  that  for  the  base-collector  junction.  The  transport  coefficients  are 

given  in  terms  of  C2lg  and  C^I^  for  the  normal  and  reverse  terms  respectively. 

Figure  13b  shows  the  graphical  interpretation  of  n and  C_I.  for  the  normal- 

© " s 

active  operation  mode.  The  total  expression  for  I is: 

B 


(b)  (60) 


The  term  (a)  in  equation  (60)  replaces  the  1/3  relationship 

CE  F 

in  Figure  2 and  the  (b)  term  replaces  I /3  . 

CC  R 

The  end  result  is  that  nine  parameters — I , V , V , C_,  I , n , C , I 

S A B * k 6 A Kr 

and  n — model  the  dc  effects  of  base-width  modulation,  recombination  and 
c 

high-level  injection  effects  in  bipolar  transistors.  Changes  in  base  transit 

time  owing  to  high-level  injection  are  incorporated  via  thr  expression: 

t. (effective)  * r Q (61) 

r f b 
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The  complete  set  of  24  Gummel-Poon  parameters  is  given  in  the  SPICE  user's 


guide  in  Appendix  III.  Of  the  24  parameters,  15  are  the  same  as  those 
needed  for  the  basic  Ebers-Moll  model.  Thus  only  an  additional  nine  must 
be  determined  to  completely  specify  the  Gummel-Poon  model. 

B.  THE  MOS  TRANSISTOR 

The  modeling  of  MOS  transistors  for  computer  implementation  is  considerably 
easier  than  is  bipolar  transistor  modeling.  The  incorporation  of  recombination 
and  high-level  injection  effects  which ^<re  necessary  for  bipolar  devices  is 
unnecessary  for  MOS  transistors.  However,  the  effect  of  bulk  charge  on 
channel  conduction  and  channel  length  (charge)  modulation  due  to  drain  voltage 
are  important  effects  which  must  be  included  in  second-order  computer  models. 
Additional  effects, which  can  be  significant, include  bias  dependent  mobility 
and  gate  capacitance.  However, these  parameters  are  taken  as  constant  in  the 
Shichman— Hodges  MOS  model  which  is  implemented  in  SPICE.  The  SPICE  model  does 
include  voltage  dependent  drain  and  source  junction  space-charge  capacitances. 

A 1/2-powwr-law  voltage  dependence  is  assumed. 


The  drain  current  for  the  MOS  ia  described  accurately  by  the  equation: 


ZuC 
K ox 

L 


2 

3 


& 


e € N . 

s ox  sub 


C 

ox 


3/2 


(62) 


where: 

G,S,D,B  - are  the  four  terminal  subscripts  and  voltages  are  positive- 
negative  as  indicated  by  the  subscripts. 
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- is  gate  oxide  capacitance  per  unit  area 


sub 

20_ 


FB 


- is  the  substrate  doping 

- is  the  strong  inversion  potential  and  is  equal  to  2kT/q  In  ^Nsub/n  j 

Qgg  ^ 

- is  the  flat-band  voltage  given  as  0 - — — where: 

MS  C 

ox 

0 - is  the  metal-semiconductor  work  function 

MS 

Q - is  the  net  surface  charge. 

SS 

The  other  terms  have  their  standard  meanings.  The  so-called  "bulk  charge" 
contribution  to  channel  conduction  is  given  by  the  second  term  in  equation  (62) 
with  its  3/2  - power  dependence.  The  gain  factor  3 can  be  modified  to  include 
channel  length  modulation  as  follows: 


ZpC 


ox 


1- 


1 

Zl7vJ 


(63) 


The  voltage  dependence  of  AL(V)  determines  the  appropriate  channel  length 
variations.  It  should  be  commented  that  mobility  variations  with  bias  can 
easily  be  Included  in  equation  (63)  by  using  a variable  p. 

The  SPICE  Shichman-Hodges  MOS  model  uses  simplified  equations  to  model 
the  general  effects  described  by  equations  (62)  and  (63).  To  model  the  "bulk 
charge"  contribution  to  conductance  the  threshold  voltage  is  modified  by  a 
voltage  dependent  term.  In  essence  the  bulk  charge  effect  is  taken  as  a constant 
during  integration  along  the  channel  so  that  a 3/2  power  dependence  is  not 
obtained.  That  is: 


- V +20  +QB 
T FB  F - — 

ox 


V +20  + v/2s  e jf 
FB  F v S OX  SUB 


* 

fv  +20  ) 

SB  ?) 


(64) 


4Y 
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If  V_-  is  defined  to  be  the  zero-V_  value  of  V„.  which  would  be  the 
TO  SB  T 

measured  threshold  for  no  "back-gate’’  bias,  then: 


V -4-Y 
TO 


(vsb*2<,f)  - (2ap)' 


(65) 


the  quantities  V , y and  20^  are  thus  input  parameters.  This  is  a mixture 
of  measured  and  calculated  parameters  since  can  be  measured  directly, 
y can  be  determined  from  a series  of  "back-gate"  bias  experiments  and  20 


is  calculated. 

The  Shichman-Hodges  model  assuaies  that  output  conductance,  G, 


DSAT 


31  /av  Iv  >V  .is  proportional  to  drain  current  in  saturation  (I 
D DS  1 DS  DSAT* 

That  is: 


DSAT 


). 


G a II 
DSAT  DSAT 


(66) 


where  lambda  is  a channel  length  modulation  parameter.  The  model  implementation 
Incorporates  the  effect  into  the'  current  expression.  If  we  define  1^  as*: 


I mSS. 

DO  2 


(VGS  “ Vt)  ” (VGD  ” Vt) 


(67) 


AF, 


the  drain  current  expression  which  corresponds  to  the  appropriate  output 
conductance  given  by  equation  (66)  is  then: 


I ■ I 


DO 


(l  * w») 


(68) 


* 

as  written, this  formulation  applies  only  above  threshold  and  for  neither 
end  of  the  channel  plnched-off.  For  normal  mode  saturation  ■ 0,  and  for 
reverse  mode  saturation  F^  m 0. 
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The  (1  + \V  ) term  given  in  equation  (68)  can  be  interpreted  as  the 

' DS ' r ^ 

4L(V) 1 


equivalent  to  the 


proportional  to  V then: 

DS 


L J 


term  in  equation  (63).  If  AL (V)  is  assumed 


- 

1 

i1  - XVDsJ 

» 1 + \V„„  for  AV  «1 
DS  DS 


(69) 


where  \ is  the  channel  length  modulation  parameter.  In  essence  the 

approximation  for  changes  in  L with  is  equivalent  to  the  linear  charge 

relationship  assumed  by  using  V for  the  Gummel-Poon  model. 

A 

In  summary,  the  MOS  model  built  into  SPICE  requires  14  parameters  for 
its  complete  specification.  The  dc  parameters  described  by  equations  (63), 
(64),  (65)  and  (66)  are  £o,  y,  V , 20  and  For  the  charge  storage 

TU  r 

elements  the  oxide— capacitances  associated  with  gate-source,  gate— drain  and 
gate-bulk  are  assumed  to  be  constant.  Junction  capacitances  of  the 
source-bulk  and  drain-bulk  have  a 1/2-power  dependence  with  user  specified 
zero-bias  values.  The  bulk  junction  potential  and  saturation  current,  as 
well  as  source  and  drain  series  resistance  are  supplied  by  the  user.  The 
complete  set  of  parameters  for  the  Shichman-Hodges  model  is  listed  in  the 
SPICE  user's  guide  in  Appendix  III. 


1 


37 


CONCLUS ION 


This  section  of  the  report  has  presented  techniques  snd 
model  formulations  which  are  basic  to  a sophisticated  computer  circuit 
analysis  program.  The  examples  have  been  taken  from  program  SPICE  which  is 
being  used  at  Stanford  for  a number  of  integrated  circuit  and  system  design 
problems  In  the  sections  which  follow,  reports  on  specific  modeling  and 
circuit  design  problems  will  be  presented.  The  emphasis  on  ac  and  transient 
analysis  capabilities  is  about  equal.  It  will  also  become  apparent  from  the 
reports  that  modeling  and  model  parameter  determination  are  major  factors  in 
determining  simulation  accuracy  and,  in  some  cases,  credibility. 
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Figure  2 Nonlinear  Bipolar  Junction  Transistor  Model  with  Charge 
Storage  Elements  and  Characteristic  Equations. 


<a> 


Figure  3.  a)  The  Exponential  Diode  Relationship  for  the  Base  Current 

due  to  I /p  with  Notation  used  for  Linearization, 
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b)  The  Linearized  Model  about  r 
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Figure  4.  The  Complete  Linearized  Bipolar  Junction  Transistor  Model 


with  Parameter  Values  to  be  Determined  as  Indicated. 


Figure  5.  The  Transistor  Circuit  from  Figure  1 for  dc 
conditions  with  Linearized  Transistor  Model. 
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Figure  6.  Voltage-controlled  Current  Source  Model  with  Reference 


Directions  and  Mode  Numbers. 
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figure  7.  A circuit  configuration  and  admittance  matrix  structure  to  illustrate 
impact  of  sparse  storage  and  reordering. 

a)  the  circuit,  uncircled  nodes  are  original  and  circled  indicate 
renumbered, 

b)  the  matrix  structure  before  reordering.  The  X's  indicate  entries, 
the  0’s  fill-in. 

c)  the  matrix  structure  after  node  reordering.  Mote  the  reduced  fill-in 


Hi) 


(tr) 


Figure  ? a)  Equivalent  Circuit  Interpretations  of  the  Integration  Formula 
for  an  Inductor  for  Forward  Euler  (F-E),  Backward  Euler  (2-E ) 


and  Trapezoidal  (TP.)  Integration  Rules. 


Figure  8 b)  Equivalent  Circuit  Interpretation  of  the  Integration  Formula, 


for  a Capacitor  for  Forward  Euler  (F-E),  Backward  Euler  (B-E) 
and  Trapezoidal  (TR) . The  square  bracketed  circuits  are  the 
equivalent  forms  appropriate  for  nodal  analysis. 
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Abstract — A survey  is  mads  of  tha  principal  techniques,  proce- 
dures. and  routines  that  are  used  in  present  programs  for  computer- 
aided  circuit  analysis.  Programs  (simulators)  are  reviewed  and 
selected  features  compared  for  the  four  maior  classes  of  circuit 
analysis:  linear  dc  and  ac.  nonlinear  dc.  nonlinear  transient,  and  linear 
pole  aero. 

I.  Introduction 

SEVERAL  recent  books  and  papers  have  cataloged 
and  compared  many  existing  and  proposed  computer- 
aided  circuit  analysis  programs  [1  J— [1 1 ).  This  paper 
is  a survey  of  the  principal  techniques  used  in  existing  pro- 
grams for  the  analysis  of  the  four  important  circuit  areas: 
linear  dc  and  ac.  nonlinear  dc,  nonlinear  transient,  and 
linear  pole  zero.  The  considerations  brought  out  are  based 
on  over  four  years  experience  in  using,  modifying,  and 
writing  more  than  twenty  programs. 

This  survey  is  organized  as  follows : first  an  overview  is 
made  comparing  nodal  and  mixed  or  state-variable  meth- 
ods of  formulating  the  circuit  equations.  Linear  dc  and  ac 
analysis  programs  are  then  considered  with  emphasis 
placed  on  techniques  for  solving  systems  of  linear  equa- 
tions. Attention  is  next  given  to  nonlinear  dc  analysis  pro- 
grams and  the  iterative  solution  of  nonlinear  algebraic 
equations.  The  extension  to  nonlinear  transient  analysis 
follows  with  an  introduction  to  several  numerical  integra- 
tion routines  used  in  the  solution  of  the  nonlinear  differen- 
tial equations.  Finally,  four  linear  pole- zero  circuit  analysis 
techniques  and  programs  are  considered  Throughout  the 
paper  a number  of  specific  references  on  computer-aided 
circuit  analysis  and  design  are  cited.  In  addition,  several 
general  computer-aided  circuit  analysis  references  are 
included  (I2J— (17], 

II.  Circuit  Equation  Formulation 
Virtually  all  presently  available  circuit  analysis  programs 
start  from  the  same  point,  an  elemental  circuit  description 
supplied  to  the  program  via  keyboard,  punched  cards,  or 
an  interactive  graphic  display  console.  This  desenption  of 
circuit  elements  and  their  interconnections  is  convened  by 
the  programs  into  a set  of  circuit  equations.  Of  interest  here 
are  the  two  major  approaches  that  are  now  used  in  formulat- 
ing these  equations.  The  first  approach  is  the  familiar  nodal 
analysis  while  the  second  is  the  mixed  or  state-variable  ap- 
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part  by  the  Joint  Services  Electronics  Program,  under  Grant  AFOSR-68- 
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Contract  DAHC04-A7-C-0031. 
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proach  which  derives  from  Bashkow's  .4-matrix  formula- 
tion (18 1. 

The  implementations  of  these  two  approaches  within 
the  present  generation  of  circuit  analysis  programs  seem  to 
have  had  a common  origin  in  the  transistor  analysis  pro- 
gram (tap)  [19 }— (201  developed  by  IBM.  tap.  though  rea- 
sonably effective,  was  short  lived  and  had  such  severe 
limitations  that  it  was  never  made  available  outside  of  I B M . 
From  tap  there  evolved  three  programs:  ecap  (21  J.  which 
used  nodal  analysis,  and  net-1  [22]  and  predict  [23], 
which  both  used  a siate-vanable  approach. 

ecap  performs  linear  dc.  linear  ac.  and  piecewise-linear 
transient  analyses  and  is  still  very  much  in  use.  It  is  con- 
sidered in  more  detail  shortly,  net-1,  a nonlinear  transient 
analysis  program,  is  available  in  the  IBM  7090/94  assembly 
languages  map  and  fap.  Therefore,  this  program  cannot  be 
used  simply  in  computer  systems  having  monitor  systems. 
net- I did.  however,  influence  the  development  of  the  pro- 
gram circus  [24].  This  program  is  considered  in  more 
detail  in  the  following  paragraphs,  predict  is  also  a non- 
linear transient  analysis  program  but  includes  no  built-in 
device  models  and  is  incapable  of  automatically  performing 
multiple  analyses  for  several  different  sets  of  element  values. 
In  addition,  it  requires  separate  analysis  runs  for  steady- 
state  fdc)  and  transient  solutions  and  suffers  from  a weak 
numerical  integration  routine.  To  overcome  these  and 
other  deficiencies,  the  sceptre  [25]  program,  which  is  con- 
sidered in  the  following  paragraphs,  was  developed. 

Nodal  Analysis  Formulation 

For  both  linear  and  nonlinear  (or  piecewise-linear)  cir- 
cuit problems,  nodal  equations  may  be  generated  by  little 
more  than  inspection.  The  relative  simplicity  of  this  ap- 
proach contrasts  sharply  with  the  manipulations  required 
by  the  state-variable  approach.  As  an  example,  consider 
the  case  of  a linear  circuit.  The  nodal  equations  are  of  the 
form 

Tr  » i 11) 

where  Y is  the  nodal  admittance  matrix,  e is  the  vector  of 
node  voltages  to  be  found,  and  i is  a vector  representing 
independent  source  currents.  The  term  y«  in  K represents 
the  sum  of  the  admittances  of  all  the  branches  connected  to 
node  i : ytl  is  the  negative  of  the  sum  of  the  admittances  of  ail 
branches  connecting  node  i and  node  j : and  it  is  the  sum  of 
all  source  currents  entering  node  k.  Thus  if  a resistor  of 
value  R connects  nodes  5 and  7.  1/R  is  added  to  v«,  and  >•-, 
and  subtracted  from  y,7  and  y7J,  while  if  a current  source  of 
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strength  l is  directed  from  node  2 to  node  3,  / is  subtracted 
from  i;  and  added  to  i9. 

Voltage  sources  are  usually  handled  in  one  of  two  ways. 
The  first  is  to  require  that  every  voltage  source  appear  in 
series  with  a resistor  so  that  the  source  may  be  transformed 
to  a Norton  equivalent  current  source.  The  second  approach 
Is  a generalization  of  the.first  but  does  not  depend  upon  a 
series  resistor.  The  approach  is  most  easily  introduced  in 
terms  of  grounded  voltage  sources.  The  nodal  equations  are 
first  assembled  including  all  elements  other  than  voltages 
sources.  Columns  of  the  admittance  matrix  corresponding  to 
grounded  voltage  source  nodes  are  then  multiplied  by  the 
value  of  the  source  and  the  result  is  subtracted  from  both 
sides  of  (1).  If  the  current  through  the  voltage  source  is 
required,  it  may  be  treated  as  an  unknown  in  place  of  the 
known  voltage.  Otherwise  the  equation  representing  the 
sum  of  the  currents  at  the  source  node  may  be  considered 
redundant  and  dropped,  thereby  reducing  the  number  of 
unknowns.  Floating  voltage  sources  and  controlled  sources 
may  be  handled  by  similar  column  and  row  operations. 

For  nonlinear  analysis,  equations  can  be  formulated  in  a 
manner  similar  to  that  used  in  the  linear  case.  The  equations 
may  be  written  symbolically  as 

»>.!)«  i<0.  (2) 

The  interpretation  of  (2)  is  that  each  equation  represents  a 
summing  of  current  contributions  at  a node.  It  should  be 
noted,  however,  that  as  now  formulated  the  nonlineanties 
are  restricted  to  be  voltage  rather  than  current  dependent. 
Fortunately,  semiconductor  devices  such  as  junction  diodes 
and  bipolar  and  field-effect  transistors  are  of  the  voltage- 
controlled  category,  and  the  above  restriction  is  usually  not 
severe. 

As  mentioned  previously,  ecap  is  based  on  nodal  analy- 
sis formulation.  Additional  nonlinear  analysis  programs 
using  nodal  analysis  include  trac  [26],  mtrac  [27], 
syscap  [28],  and  bias-3  [29],  while  linear  programs  based 
on  a nodal  approach  include  natfmqs  [30],  [31  ],  and  from 
lisa  [32],  [33]  the  programs  acca,  poly,  and  trfn. 

State-  Variable  Analysis  Formulation 

The  primary  reason  for  using  the  state-variable  approach 
is  that  it  yields  a set  of  first-order  linear  or  nonlinear  differ- 
ential equations  in  a minimal  set  of  unknowns.  Deriving 
a set  of  such  equations  explicitly  seemed  a desirable  objec- 
tive for  use  with  earlier  numerical  integration  algorithms. 
As  is  brought  out  in  the  following  paragraph,  this  is  no 
longer  essential. 

From  an  elementary  viewpoint,  the  state-variable  for- 
mulation proceeds  as  follows  [34]:  a proper  tree  consisting 
of  all  voltage  sources,  as  many  capacitive  branches  and  as 
few  inductive  branches  as  possible,  and  no  current  sources 
is  selected.  The  state  variables  are  chosen  to  be  the  capacitive 
tree-branch  voltages  or  charges  and  inductive  tree-link  cur- 
rents or  fluxes.  A fundamental  cut-set  equation  is  con- 
structed for  each  capacitive  tree  branch  and  a fundamental 


loop  equation  for  each  inductive  tree  link.  Independent 
sources  are  withdrawn  as  separate  terms  in  the  equations 
and  the  equations  are  normalized  with  respect  to  the  ele- 
ment values  associated  with  each  state  variable. 

For  a linear  circuit,  the  preceding  procedure  results  in  a 
matrix  equation  of  the  form  1351 

x = Ax  +•  8u  (3) 

where  A is  a coefficient  matrix  relating  the  state  vector  x 
to  its  derivative  x and  B is  a coefficient  matrix  coupling  the 
effects  of  the  independent  source  vector  a.  Any  other  de- 
sired network  variables  can  be  expressed  in  terms  of  the 
state  variables  and  independent  sources.  For  the  linear 
case,  this  procedure  results  in  an  equation  of  the  form 

y » Cx  + Du  (4) 

where  y is  a vector  of  desired  output  vs  nables  and  C and  D 
are  again  coefficient  matrices. 

For  the  nonlinear  case.  ( 3 1 and  (41  are  of  the  form 

x =•  A(x.  u.  r)  (5) 

y = C(x,  a,  t).  16) 

An  alternative  description  is  based  on  extracting  (he  non- 
linear portion  of  the  circuit  from  the  linear  portion.  This 
allows  the  state  equations  1 5)  and  (6)  to  be  written  in  the 
following  equivalent  form  (15  ] : 


Ax  + Bu  + B’w 

(7) 

Cx  + Du  + D’w 

(8) 

Ex  + Fu  +■  F'». 

|9) 

Again,  A,  B,  B',  C.  D,  D‘,  E.  F,  and  F'  are  coefficient  matrices, 
while  x and  ■ remain  the  state  variables  and  independent 
source  vector,  respectively.  The  vector  w represents  the  con- 
trolling or  independent  variable  set  associated  with  the 
noniinearitics  gfw)i  In  the  state-vanable  approach,  no 
assumption  is  made  regarding  voltage  or  current  control. 

While  the  equation  formulation  procedure  outlined  above 
applies  in  general,  details  vary  from  program  to  program. 
Both  emeus  and  sceptre  emp'.ov  variations  of  Bryant’s 
method  [36],  [37]  as  modified  by  Wilson  and  Massena  [38  ]. 
Other  programs  relying  on  the  state-variable  approach  in- 
clude cornap  [39],  [40],  belac  [4 1 ].  and  cirpac  [42],  (43 1. 
In  all  cases,  the  formulation  procedures  entail  extensive 
matrix  manipulations.  The  bookkeeping  associated  with 
these  manipulations  usually  requires  considerable  amounts 
of  core  storage.  Loops  of  capacitors  and  voltage  sources 
and  cut-sets  of  inductors  and  current  sources  are  accommo- 
dated by  additional  manipulations  of  the  state  equations. 

III.  Linear  DC  and  AC  Analysis 

The  first  of  the  four  circuit  analysis  areas  to  be  treated  is 
linear  dc  and  ac  analysis.  The  desired  outputs  are  the  various 
voltages  and  currents  within  a circuit,  possibly  as  a function 
Of  frequency.  As  virtually  all  programs  of  this  category  use 
a nodal  analysis  formulation,  it  is  assumed  that  the  circuit 
equations  to  be  solved  have  the  form  of  ( I ).  The  major  tech- 


16 


IEEE  TRANSACTIONS  ON  ClMCV IT  THEORY,  JANUARY  197 1 


niques  presently  used  to  solve  such  systems  of  linear  equa- 
tions include  Gaussian  elimination,  pivoting,  LU  factoriza- 
tion. and  sparse  matrix  methods.  These  techniques  are 
outlined  in  the  following  paragraphs  and  are  followed  by 
a comparison  of  the  features  of  several  available  analysis 
programs. 

A dc  analysis  can  be  considered  the  special  case  of  an  ac 
analysis  performed  at  zero  frequency.  However,  such  an 
approach  is  usually  not  used  for  two  reasons:  first,  if  the 
formulation  proceeds  on  a nodal  admittance  basis,  the 
infinite  admittance  represented  by  an  inductor  at  zero  fre- 
quency requires  special  consideration.  Second,  where  only 
resistive  elements  and  independent  and  controlled  sources 
are  present,  the  circuit  equations  contain  only  real  coeffi- 
cients and  hence  may  bf  olved  on  a computer  using  only 
real  variables.  If  complex  variables  are  used,  the  computa- 
tion time  may  be  increased  by  as  much  as  a factor  of  4 be- 
cause multiplication  and  division  (both  considered  long 
operations  as  compared  with  the  short  operations  of  addi- 
tion and  subtraction)  require  more  time  than  for  long 
operations  with  real  variables. 

Guussian  Elimination 

With  the  distinction  in  mind  that  the  system  of  nodal 
equations  (1)  involves  only  real  variables  in  the  dc  case  and 
complex  variables  in  the  ac  case,  the  solution  of  (11  can  be 
written 

r «*  Y~'i  (10) 

where  ¥~ 1 is  thr  inverse  of  the  nodal  admittance  matrix. 
Computationally,  the  efficiency  of  this  approach  can  be 
evaluated  in  terms  of  the  number  of  long  operations  re- 
quired. It  can  be  shown  that  the  number  of  long  operations 
required  to  invert  an  n x n matrix  is  it3,  while  the  number  of 
long  operations  required  to  compute  its  oroduct  with  a vec- 
tor of  dimension  n is  n2 . Thus  the  total  number  of  long  oper- 
ations required  by  this  approach  is 


The  number  of  long  operations  required  to  obtain  a solu- 
tion to  1 1 ) can  be  reduced  by  a factor  of  3 using  Gaussian 
elimination  (44). 

This  procedure  consists  of  two  steps.  The  first  step  con- 
sists of  converting  the  nodal  admittance  matnx  Y to  an 
equivalent  upper  triangular  matrix,  i.e..  a matrix  with  all 
zero  elements  below  the  diagonal.  The  second  step  which  is 
referred  to  as  back  substitution,  consists  of  solving  the  nth 
equation  containing  only  r.  for  t>„,  the  n-  1st  equation  for 
e,. , in  terms  of  c„  etc.  Gaussian  elimination  is  illustrated 
in  1 1 1 ) for  a third-order  system 
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from  S'?  and  I AH/ Ail)  S'?  from  S'?.  Symbolically  this 
transformation  can  be  represented  by  the  equations 
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which  yields  the  system 
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Finally,  the  unknown  e2  is  eliminated  from  S'?  by  sub- 
tracting IlMV:11  from  thus  obtaining 

jfl3l  _ mlt  API  1 1 

(5  j s 0 | — “ [ 

s'?  - r? 

s'?  - s'?  - j|U2'  (14) 

resulting  in  the  triangularized  system 


(15) 


This  completes  the  first  step  in  the  Gaussian  elimination 
procedure.  Back  substitution  is  now  performed  to  obtain 
the  final  solution  as  follows : 
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(III 


where  the  superscript  I indicates  the  initial  system  of  equa- 
tions. As  indicated  previously,  the  unknown  r,  is  eliminated 
from  equations  f?  and  S'?  by  subtracting  (All!  Atlas'? 


Careful  enumeration  shows  that  for  an  nth  order  system 
the  number  of  long  operations  required  by  Gaussian  elimi- 
nation is 

n 3 , n 

7 + " -3' 

LU  Transformation 

A modification  of  Gaussian  elimination  which  is  useful 
when  more  than  one  source  or  right-hand-side  vector  i is 
to  be  considered  is  the  LU  transformation  (44).  This  pro- 
cedure consists  of  partitioning  Y into  an  upper  triangular 
matrix  U and  a lower  triangular  matrix  L (usually  with  ones 
on  the  diagonal)  such  that 

LU  =«  Y.  (17) 

This  technique  is  illustrated  shortly.  The  resulting  system  is 
solved  in  two  stages.  First. 

Uv=lT'  i = i*  (1*1 
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and  secondly. 


v 


ir  li*. 


(19) 


In  this  case  since  both  L and  U are  triangular,  IT ' and  U~ ' 
are  trivial  to  compute.  Note  that  ( 19)  is  the  same  back  sub- 
stitution performed  during  Gaussian  elimination,  while  (18) 
is  also  equivalent  to  back  substitution.  To  obtain  V the 
same  reduction  procedure  is  performed  as  in  Gaussian 
elimination.  Thus,  in  terms  of  the  earlier  third-order 
example. 


A'i  A'i 
y'ii  AH 
0 yft 


(20) 


Similarly,  L is  found  to  be  a byproduct  of  the  triangulariza- 
tion  procedure.  For  the  example. 
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Ail 


AH 

Wi 


0 

0 

1 


(21) 


formed  in  double  precision  at  the  expense  of  only  one  ad- 
ditional double  precision  vector  of  dimension  n.  The  only 
major  added  expense  is  computation  time. 

The  most  commonly  used  approach  of  reducing  error 
(used  by  Ralston  with  the  preceding  method)  is  pivoting  [4$]. 
Partial  pivoting  amounts  to  scanning  the  elements  of  the  ith 
column  of  (he  matrix  below  the  diagonal  at  the  it h step 
tn  the  reduction  and  determining  the  element  of  largest 
magnitude.  The  row  of  this  element  is  exchanged  with  the 
ith  row  and  the  reduction  is  continued.  For  the  example 
considered  previously,  if  at  the  second  step  ly&’j  > |/g|,  y'}i 
would  be  chosen  as  the  pivot  element,  S'}'  and  S'}'  would  be 
exchanged  and  y'ii  would  be  eliminated  from  S'}'  to  obtain 
S'}'.  This  technique  tends  to  reduce  the  magnitude  of  the 
term  being  subtracted  at  each  step  and  thus  improve  ac- 
curacy. Further,  it  preserves  column  order  and  hence  the 
order  in  the  solution  vector  v. 

Complete  pivoting  involves  finding  the  largest  element  in 
the  as  yet  unreduced  (ft  — i)  x (n  — i)  submatrix  at  the  ith  step, 
exchanging  rows  and  columns  such  that  it  appears  on  the 
diagonal  in  the  ith  row  and  column.  Since  column  order  is 
not  preserved  except  at  the  expense  of  additional  book- 
keeping, complete  pivoting  is  seldom  used. 


The  elements  of  L are  the  coefficients  in  (12)  and  (14). 

A convenient  aspect  of  the  preceding  procedure  is  that 
since  the  diagonal  elements  of  L are  known  to  be  l's,  L and 
V may  share  the  same  memory  locations  originally  assigned 
to  F.  Note  further  that  the  long  operations  count  is  the  same 
as  for  Gaussian  elimination.  Once  L and  U have  been  com- 
puted, (18)  and  (19)  can  be  applied  repeatedly  for  different 
source  vectors  1.  For  m source  vectors,  the  total  long  opera- 
tions count  becomes 


Pivoting 

Further  methods  of  reducing  core  requirements  and 
operations  counts  are  brought  out  below.  First,  however, 
mention  should  be  made  of  error  control.  As  can  be  seen 
from  1 12)  and  ( 14),  a division  and  a subtraction  are  required 
at  each  step  in  the  triangulanzation  process.  Suppose  a 
computer  which  represents  a number  accurately  to  d digits, 
performs  arithmetic  operations  correctly  to  2d  digits,  and 
then  rounds  the  result  to  U digits.  It  can  be  shown  [44]  that 
the  absolute  error  resulting  from  such  arithmetic  operations 
is  given  by  l+<h  10"'  where OtSipsS.  In  the  floating-point 
notation,  the  operations  of  subtraction  and  division  yield 
results  where  magnitudes  are  less  than  those  of  the  operands. 
The  relative  error  in  the  result  is  larger;  thus  division  and 
subtraction  steps  reduce  accuracy.  The  lower  accuracy  is 
felt  more  severely  on  computers  with  smaller  word  lengths 
where  d is  reduced.  One  obvious  solution  is  to  declare  all 
variables  to  be  double  precision  and  perform  all  operations 
in  double  precision.  However,  this  significantly  increases 
both  time  and  core  memory  requirements. 

One  compromise  which  can  be  made  is  described  by 
Ralston  145  J.  He  shows  that  the  critical  steps  involving 
divisions  and  accumulation  of  partial  sums  can  be  per- 


Sparse  Matrix  T echniques 

In  view  of  the  n3  dependence  of  the  long  operations  count, 
computation  time  can  be  expected  to  increase  significantly 
as  larger  circuits  with  more  nodes  are  considered.  Several 
recent  papers  (46)- [48]  have  focused  attention  on  taking 
advantage  of  sparsity  in  the  nodal  admittance  matrix.  There 
are  basically  three  associated  economies.  First,  efficient 
means  have  been  found  by  which  only  the  nonzero  entries 
of  the  matrix  need  be  stored,  thus  effecting  a savings  in  core 
memory.  Second,  it  is  possible  to  process  only  the  nonzero 
entries  at  each  step  in  the  triangular  reduction.  Finally,  the 
order  in  which  variables  are  eliminated  can  be  chosen  to 
preserve  sparsity.  The  long  operations  count  and  computa- 
tion time  are  (hen  reduced.  This  savings  becomes  even  more 
significant  when  the  same  equations  must  be  solved  many 
times,  as  in  multifrequency  analysis.  The  optimal  order  for 
eliminating  variables  need  only  be  determined  once. 

By  way  of  illustration  (86),  in  a new  program  developed 
at  the  University  of  California,  Berkeley,  by  Prof.  R.  Rohrer 
and  his  students,  the  third  technique  of  optimal  ordering 
together  with  nonzero  storage  leads  to  a long  operation 
count  more  closely  proportional  to  n than  to  n3.  In  the 
analysis  of  a typical  operational  amplifier  of  22  nodes,  the 
number  of  long  operations  was  reduced  from  2660.  using 
Gaussian  elimination,  to  130. 

Linear  DC  and  AC  Analysis  Programs 

Many  linear  analysis  programs  are  available  including 
GCap.  the  acca  portion  of  lisa,  and  rohrerx.'  As  previ- 
ously mentioned,  these  programs  are  based  on  a nodal 
analysis  formulation.  Both  acca  and  rohrerx  have  free 

1 rohrerx  is  a proeram  written,  xml  developed  in  IVOR  by  the  1C 
Group.  Electronics  Research  Laboratory.  University  of  California. 

Berkeley. 
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format  input  languages  as  do  certain  time-sharing  versions 
of  ecap.  The  primary  advantages  of  ecap  are  that  it  is  well 
documented  [SO  | and  is  available  through  many  commercial 
time-sharing  companies.  Us  disadvantages  are  that  it  is 
cumbersome  to  modify  and  the  original  input  language  is 
cumbersome  to  use. 

acca,  as  part  of  the  larger  lisa  package,  shares  a general 
input  routine  and  has  the  advantage  of  being  able  to  inter- 
act somewhat  with  other  portions  of  lisa.  Its  disadvantage 
is  that,  if  used  as  part  of  lisa,  it  is  expensive  to  load. 
ROHRHtx  does  make  partial  use  of  sparse  matrix  techniques 
and  is  therefore  relatively  fast.  It  presently  suffers,  how- 
ever, from  a lack  of  documentation. 

At  the  present  time  efforts  are  being  made  to  decrease 
computation  time  and  to  incorporate  sensitivity  and  noise 
performance  using  the  adjoint  network  [51  ].  This  work  is 
of  major  importance  in  the  development  of  several  com- 
puter-aided circuit  design  packages  (52]-[54). 

IV.  Nonuneaji  DC  Analysis 

Most  nonlinear  dc  analysis  programs  are  incorporated 
into  more  general  transient  analysis  programs.  Equation 
formulation  approaches  are  divided  between  nodal  and 
state-variable  analysis.  In  either  approach,  the  nonlinear  dc 
analysis  problem  reduces  to  one  of  solving  simultaneous 
nonlinear  algebraic  equations.  As  described  in  the  following 
paragraphs,  iterative  methods  are  used.  Of  particular  im- 
portance is  the  convergence  properties  of  the  solution  algo- 
rithm. Several  approaches  which  are  used  in  different  pro- 
grams to  improve  convergence  are  examined  in  this 
section.  Finally,  various  programs  are  again  compared. 

Functional  Iteration 

The  solution  method  used  by  virtually  all  of  the  presently 
available  nonlinear  dc  analysis  programs  is  based  on  the 
Newton- Raphson  iteration  technique.  It  is  one  of  a broad 
class  of  techniques  known  collectively  as  functional  itera- 
tion methods  [44). 

Given  a set  of  nonlinear  equations  of  the  form  of  (2),  for 
the  dc  case  the  equations  may  be  expressed 

g(sl  - 1.  (22) 

The  solution  technique  is  to  start  from  some  initial  set  of 
values  and  to  generate  a sequence  of  iterates  - - • ", 

s'".  «"♦",••  ■ which  converge  to  the  solution  ». 

Newton-Raphson  iteration  is  most  easily  introduced  by 
considering  the  case  of  a single  nonlinear  equation 

«<v)  - 0.  (23) 

The  function  g(r)  can  be  expanded  about  some  point  o0  in  a 
Taylor  series  to  obtain 

g<®)  - g( »„)  + (o  - o0)g'(o0)  + • ■ • » 0 (24) 

where  the  prime  denotes  differentiation  with  respect  to  o. 

If  only  first -order  terms  are  retained,  a rearrangement  of 
(24)  yields 


o 


Vo  ~ 


g0>  o) 
9'(v  o) 


(25) 


The  form  of  (25)  suggests  that  a sequence  of  iterates  might  be 
generated  by  the  following: 


Equation  (26)  is  the  Newton-Raphson  iteration  function 
for  the  scalar  case.  Note  that  at  the  solution  5.  9lC)-0  and 
!)>’*  " = (/"  as  would  be  expected.  The  geometrical  interpre- 
tation of  (26)  is  illustrated  in  Fig.  I for  the  simple  case  of  a 
current  source  driving  an  ideal  diode.  The  line  tangent  to 
the  nonlinearity  at  the  point  (iM,  giv'"))  has  the  slope 
g'(vt’,X  Its  intercept  with  the  voltage  axis  defines  the  next 
voltage  iterate  in  the  sequence  as  shown  in  the  figure. 

The  generalization  of  the  Newton-Raphson  procedure 
to  a system  of  n equations  is  given  by 


. ,>«  _ J-VW*)  (27) 


where  the  Jacobian  Jit)  of  the  function  g(r)  is  given  by 


dg  i dg,  dg, 

dp,  dv2  do , 

•/(»)  = 

dg. dg. 

do,  civ 


(28) 


A physical  interpretation  of  the  elements  of  J(t)  is  brought 
out  below. 

The  direct  application  of  (27)  necessitates  computing  the 
inverse  of  the  n x n Jacobian  matrix.  As  indicated  previously, 
the  operation  count  for  inverting  a matrix  and  multiplying 
the  result  by  a vector  is  n2  + n2. 

An  alternative  procedure  for  obtaining  new  iterates  is  to 
solve  the  linear  system  of  equations 

J(s>")(»1*'  - t1**")  = jl*1”).  (291 

A second  alternative  procedure  often  used  with  nodal 
analysis  is  to  employ  the  system  of  equations 

J(W  ii  _ J(fit)#l«  _ (30) 

The  right-hand  side  of  (30)  is  found  to  have  a particularly 
simple  interpretation  as  also  brought  out  below. 

Gaussian  elimination  applied  to  either  (29)  or  (30)  re- 
duces the  long  operation  count  to  (n3/3)  + nJ— |n/3).  An 
even  greater  advantage  is  obtainable  using  sparse  matrix 
methods.  The  locations  of  the  nonzero  elements  of  the 
Jacobian  matrix  are  fixed  by  the  circuit  topology  and  remain 
unchanged  from  iteration  to  iteration.  The  additional  time 
required  on  the  first  iteration  to  record  the  nonzero  struc- 
ture and  determine  the  optimal  variable  elimination  order 
is  small  in  comparison  to  the  total  computation  time  re- 
quired as  the  number  of  iterations  becomes  large. 

Newton- Raphson  Applied  To  Nonlinear  DC  Analysis 

A physical  interpretation  of  the  Jacobian  matrix  and  (he 
Newton-Raphson  method  can  be  made  using  the  diode 
circuit  of  Fig.  2(a).  The  exponential  nonlinearity  of  the 
diode  is  linearized  about  some  trial  solution  voltage  Tn. 
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Fig.  1.  Ncwlon- Raptoon  itenuon. 


l(j)  V 


Fig.  2.  Nonlinear  dc  anaiyu*.  (a)  Diode  circuit,  (b) 

approximation.  (c)  Unearned  circuit  model. 


Unearned  diode 


This  is  equivalent  to  a Taylor  series  expansion  as  indicated 
previously  where  only  first-order  terms  are  retained.  The 
expansion  is  of  the  form 


lo 


Id  + («'- 


lv»v„ 


= l,(er«rT  - I)  + (V  - Vt)h-t 

yT 

m I oo  + ( V ~ K>)0oo 


Voir-r 


(31) 

(32) 

(33) 


where  /„ „ is  recognized  as  the  current  through  the  diode 
lurtvspondmg  to  the  voltage  K>,  and  gM  is  recognized  as 
the  dynamic  conductance  corresponding  to  the  voltage 
V0.  Since  the  diode  characteristic  as  described  by  (33) 
has  now  been  linearized,  the  diode  may  be  modeled  in  terms 
of  a Norton  equivalent  current  source  l0oo  111  parallel 
with  the  conductance  goo-  As  can  be  seen  from  Fig.  2(b), 
/two  is  given  by 

loo o m Ido  ~ 9o o^o-  (34) 

Hence;  1 32)  may  be  written  in  terms  of  / D*0  as 

lom  9 do  y + It mo-  (35) 

The  nodal  equation  for  the  complete  linearized  circuit  of 
Fig.  21c) is 

(0  + 9„o)y  - I - 'd oo-  (36) 


In  terms  of  iterate  values 

(G  + 9'J')F<”'>  =./-/&.  (37) 

If  (30)  is  compared  with  (37)  the  physical  interpretations 
of  the  Jacobian  and  the  right-hand  side  of  (30)  become  more 
apparent.  The  Jacobian  consists  of  the  nodal  conductance 
matrix  of  the  linear  elements  of  the  circuit  together  wirh  the 
linearized  conductances  associated  with  each  nonlinear  cir- 
cuit element.  The  vector  on  the  right-hand  side  of  (30) 
consists  of  independent  source  currents  and  the  Norton 
equivalent  source  currents  associated  with  each  nonlinear 
circuit  element.  Thus  at  each  iteration  in  a nodal  analysis, 
the  linearized  conductances  and  Norton  equivalent  source 
currents  must  be  recomputed  and  the  linearized  nodal 
conductance  equations  reassembled. 

The  application  of  Newton-Raphson  iteration  lo  the 
state-variable  approach  is  also  straightforward  once  all 
required  coefficient  matrices  have  been  assembled.  Recall 
that  the  equations  formulated  by  the  state-variable  method 
can  be  put  in  the  form  of  (7H9).  The  state  vector  x may 
represent  capacitor  voltages  and  inductor  currents  in  which 
case  the  steady-state  dc  solution  is  characterized  by  i »o. 
This  corresponds  lo  setting  the  current  through  capacitors 
and  voltages  across  inductors  to  zero.  Next,  (7)  may  be 
solved  for  x in  terms  of  » and  the  result  substituted  into  |9) 
to  obtain  a system  of  nonlinear  equations  in  terms  of  w. 
Once  w has  been  obtained  x and  y can  also  be  obtained. 
Note  that  since  the  coefficient  matrices  involved  are  con- 
stant. the  initial  manipulation  necessary  lo  obtain  the  single 
system  of  nonlinear  equations  in  * need  only  be  performed 
once. 

Comergence 

For  the  nonlinear  dc  analysis  approach  just  outlined,  the 
problem  of  convergence  (o  a solution  is  now  considered. 
Proofs  that  an  algorithm  will  converge  depend  upon  a 
priori  knowledge  of  an  initial  guess  sufficiently  close  to  the 
solution.  Because  this  knowledge  is  usually  noi  present  in  a 
nonlinear  dc  analysis,  all  techniques  incorporated  into 
analysis  programs  for  improving  convergence  of  the 
Newton-Raphson  iteration  technique  are  supported  by 
purely  empirical  justification. 

The  exponential  nonlineanties  usually  associated  with 
diodes  and  bipolar  transistors  are  single-valued  mono- 
tomcally  increasing  continuous  functions.  However,  these 
expressions  are  strongly  increasing  functions.  For  large 
reverse  bias,  the  slope  approaches  zero,  while  for  large 
forward  bias  the  exponential  tends  to  infinity.  Convergence 
may  be  slowed  or  the  iteration  procedure  may  be  stopped 
when  numbers  exceed  a computer-imposed  limit.  An  ap- 
proach commonly  used  to  prevent  such  overflows  is  illus- 
trated in  Fig.  3 again  for  a simple  diode  characteristic.  As 
previously  indicated,  at  a trial  junction  voltage  if”,  the 
characteristic  is  modeled  in  terms  of  a linearized  approxi- 
mation as  shown.  A new  iterate  value  if" * " greater  than 
if”  must  correspond  to  a solution  on  the  linearized  char- 
acteristic at  the  point  I'.  Two  choices  for  a new  trial  operat- 
ing point  are  immediately  available.  The  first  is  to  update 
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Fif.  3.  Selection  of  new  trt*l  operating  point. 

with  voltage  by  proceeding  vertically  to  a new  point  on  the 
exponential  characteristic.  This  is  the  standard  Newton- 
Raphson  iteration  and  can  be  carried  out  by  a straight- 
forward evaluation  of  the  exponential  function.  The  second 
choice,  which  prevents  the  overflow  problem,  is  to  update 
with  current  by  moving  back  horizontally  to  the  exponential 
characteristic.  This  is  done  by  computing  (he  logarithm  of 
the  current  corresponding  to  the  point  I'  and  results  in  the 
selection  of  point  I as  (he  new  trial  linearization  point.  Note 
that  this  procedure,  which  is  particularly  suited  for  junc- 
tion diode  and  bipolar  transistor  circuits,  can  be  used  with 
both  nodal  analysis  and  state-variable  analysis.  This  modi- 
fied Newton-Raphson  method  can  be  used  with  other  non- 
linear devices  but  does  require  that  they  are  described  by 
functions  which  have  an  explicit  inverse  and  for  which  both 
the  function  and  its  inverse  are  single-valued  and  con- 
tinuous. For  diodes  and  transistors  the  vanishing  slope  in 
the  reverse  direction  is  usually  handled  by  placing  a small 
leakage  conductance  across  each  junction. 

Broyden  [55]  has  recently  proposed  a variation  of  the 
Newton-Raphson  technique.  The  method  incorporates 
two  modifications.  The  norm  of  the  vector  giek  which  must 
tend  to  zero  as  a solution  is  approached,  is  never  allowed  to 
increase.  The  method  also  avoids  computation  of  the  in- 
verse Jacobian  matrix  at  each  iteration.  Instead,  an  arbitrary 
approximation  to  the  matrix  is  chosen  at  the  first  iteration 
and  then  successively  updated.  Branin  and  Wang  [36]  have 
applied  the  method  to  nonlinear  dc  problems,  including 
statistical  analysis.  More  recently.  Broyden  [57]  has  con- 
cluded that  enforced  norm  reduction  is  not  always  ad- 
vantageous and  that  an  adaption  of  his  method  used  in 
conjunction  with  a particular  form  of  Davidenko's  method 
(58 ] may  converge  more  rapidly.  In  addition.  Brown  [59] 
has  proposed  a variation  of  the  Newton-Raphson  technique 
in  which  an  inequality  constraint  is  applied  to  the  norm  of 
the  vector  -•/“'(»)$(»). 

Three  termination  criteria  are  commonly  used.  The  first 
is  to  stop  iterating  when  the  absolute  difference  between 
each  unknown  voltage  or  current  iterate  and  its  previous 
value  is  reduced  below  some  preset  minimum.  The  second 
is  to  stop  when  the  relative  error,  defined  as  the  absolute 
difference  divided  by  the  value  of  the  iterate,  is  reduced  to  a 
preset  minimum.  Finally,  an  approach  sometimes  used  with 
nodal  analysis  is  to  require  that  the  sum  of  the  currents  at 
each  node  be  reduced  to  a preset  minimum.  All  of  the  ap- 
proaches have  advantages  and  disadvantages  in  specific 
cases  and  none  is  superior  in  general. 


.Won linear  DC  Analysis  Programs 

Both  net- I and  its  predecessor  tap  formulate  the  dc 
equations  separately  from  the  transient  situation.  Further, 
both  programs  consider  the  linear  and  nonlinear  portions 
of  a network  separately  for  dc  analysis.  The  linear  portion 
is  analyzed  using  a modification  of  Kron's  method  of 
tearing  [60  ].  Nonlinearities  are  treated  as  a set  of  side  con- 
straints. and  Newton-Raphson  iteration  is  used. 

In  sceptre  different  trees  are  chosen  for  dc  and  transient 
analyses.  A modified  Newton-Raphson  procedure  [61  ] is 
used  which  ensures  that  the  true  operating  point  on  an 
exponential  is  approached  from  below.  This  amounts  to 
updating  with  current  whenever  junction  voltage  increases. 

Both  circus,  using  the  state-variable  approach,  and 
TRAC,  using  nodal  analysis,  employ  Newton-Raphson 
iteration  while  updating  with  voltage  in  the  third  quadrant 
and  updating  with  current  in  the  first  quadrant,  bias-3. 
using  nodal  analysis  and  Newton-Raphson  iteration,  up- 
dates with  current  in  the  first  quadrant  when  a junction 
voltage  increases  and  with  voltage  when  it  decreases.  In  the 
third  quadrant  bias-3  always  updates  with  voltage  while 
modeling  the  exponential  characteristic  in  terms  of  a line 
through  the  origin  rather  than  a tangent.  This  does  not 
affect  the  dc  solution  and  yet  insures  that  if  a junction  volt- 
age becomes  positive,  the  new  trial  linearization  point  will 
lie  in  the  first  quadrant. 

trac.  circus,  and  bias-3  converge  reasonably  well  on 
most  moderately  sized  circuits  of  up  to  40  nodes.  Nonlinear 
models  built  into  trac  include  junction  diodes  and  bipolar 
transistors.  The  program's  input  format  is  cumbersome 
while,  as  supplied  by  the  Harry  Diamond  Laboratory. 
trac  includes  several  assembly  language  subroutines  which 
may  require  translation,  circus  has  the  advantage  of  a free 
format  input  language,  zener-diode,  tunnel-diode,  unijunc- 
tion. and  junction  field-effect  transistor  models  in  a stored 
library.  It  is.  however,  a much  more  complex  program  to 
implement  and  modify,  sias-3  is  relatively  small  but  is 
limited  to  nonlinear  dc  analysis  of  bipolar  transistor  circuits. 

sceptre  is  by  far  the  most  flexible  program  in  that  models 
may  be  built  by  the  user,  nested,  and  recalled  as  necessary. 
The  price  of  this  flexibility  is  size  and  complexity  as  the  pro- 
gram consists  of  over  15  000  statements.  Its  dc  solution  has 
also  been  known  to  suffer  from  poor  convergence  proper- 
ties. 

One  recent  program  should  also  be  mentioned.  The 
dicap  portion  of  syscap  is  an  extremely  general  dc  analysis 
program.  It  is  large  and  is  currently  available  only  to  users 
of  Control  Data  Corporation's  Cybernet  remote  batch 
system. 

V.  Nonlinear  Transient  Analysis 

The  general  procedure  for  the  transient  analysis  of  a > 
nonlinear  circuit  is  to  evaluate  the  state  of  the  circuit  at  a 
given  point  in  time  and  to  extrapolate  ahead  to  a new  time, 
point. 

The  computation  time  required  for  such  an  analysis  pro- 
gram is  directly  proportional  to  the  number  of  time  incre- 
ments into  which  the  analysis  (simulation!  lime  must  he 
divided.  The  total  simulation  time  is  usually  a multiple  of 
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the  largest  tune  constant  {smallest  eigenvalue!  associated 
with  the  circuit  linearized  at  a time  point.  On  the  other 
hand,  for  a large  class  of  programs  Including  net- I,  circus. 
and  sceptre,  the  length  of  a time  step  A in  order  to  retain 
solution  accuracy  is  determined  by  the  smallest  time  con- 
stant (largest  eigenvalue)  within  the  circuit.  Physically  in- 
terpreted, this  implies  that  if  high-frequency  devices  are 
used  in  a low-frequency  application,  the  computation  time 
for  an  adequate  simulation  may  be  excessive.  To  understand 
how  this  problem  is  alleviated  in  programs  such  as  traC 
and  cirpac.  it  is  necessary  to  examine  the  numerical  inte- 
gration process  and  the  associated  problem  of  stability. 

Before  proceeding,  however,  the  use  of  the  term  numerical 
integration  should  be  clarified.  Strictly  speaking,  integration 
is  associated  with  finding  the  area  under  some  known  func- 
tion. The  nonlinear  transient  analysis  problem  is  one  of  ob- 
taining the  solution  to  a set  of  first-order  nonlinear  differ- 
ential equations  of  the  form  (S).  For  the  nonlinear  situation. 
(5)  can  be  written 


The  notation  indicating  the  dependence  of  /(x(t))  on  uU) 
has  been  dropped  for  convenience.  Some  initial  condition 


must  be  specified  a priori  or  found  as  the  result  of  a dc 
analysis.  Because  of  the  similarity  of  the  formulas  used  in 
obtaining  a solution  to  (38)  to  numerical  integration  for- 
mulas, the  term  has  come  to  be  looseiy  applied  to  both 
numerical  processes. 

Explicit  Integration 

Suppose  at  some  point  in  time  that  x(r)  is  approximated  by 
+ (40) 

If  (41)  is  substituted  into  (38),  the  result  can  be  rearranged 
into  the  form 

xlt  + A)  = A/(x(r))  + x(r).  (41) 

Here  xlt+h)  is  defined  explicitly  in  terms  of  x[t).  The  nu- 
merical integration  procedure  defined  by  (41)  is  known  as 
the  forward  Euler  method  (16]  and  is  illustrated  in  Fig.  4(a). 
If  i(r)  is  the  exact  solution  to  (38),  it  is  easy  to  visualize  the 
gross  errors  which  can  result  from  choosing  A too  large. 
The  similarity  of  the  right-hand  side  of  (41)  to  a truncated 
Taylor  series  expansion  about  t suggests  that  the  error  in 
x(r  + A)  will  be  proportional  to  h2  multiplied  by  the  second 
derivative  of  x|r)  evaluated  someplace  between  t and  t+h. 

The  forward  Euler  method  defined  by  (41)  is  so  low  in 
accuracy  that  it  is  seldom  used.  Nonetheless,  it  serves  to 
illustrate  the  idea  that  while  the  exact  solution  to  (38)  must 
be  continuous,  the  computed  solution  is  at  best  a piecewise- 
linear  approximation  to  the  exact  solution. 

In  Fig.  44a)  the  difference  between  the  exact  solution  and 
the  computed  solution  suggests  the  plausibility  of  using  the 
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Fig.  4.  Numerical  integration,  (a)  Forward  Euler  method. 

(b>  Improved  predictor-corrector  method. 

average  value  of  the  derivatives  at  i and  t + h in  (41).  This 
results  in  what  is  called  the  trapezoidal  integration  method 

x(r  + ft)  = x(r)  + * [/(x(r))  + f(Mt  + A))].  (42) 

Equation  (41)  can  be  used  to  provide  an  initial  approxim-.- 
tion  to  x(t+h)  and  thus  /(xlt+h)).  This  is  made  cl r,c,  if 
(41)  and  (42)  are  rewritten  as  follows : 

x,(t  + h)  =•  A/(xt(t))  + xc(t)  (43) 

x,(i  + h)  =■  xc(f)  + ^ [f(xr!'!)  + /(x,(t  + A))].  (44) 

Equations  (43)  and  (44)  constitute  a predictor-corrector 
pair,  the  predictor  equation  (43)  providing  xr(t  + h ) as  an 
explicit  function  of  x^t)  and  the  corrector  equation  (44) 
providing  xc(t+h)  as  an  explicit  function  of  xc(r)  and 
x,(r  + A).  This  method  is  illustrated  in  Fig.  4(b).  The  increased 
accuracy  which  can  result  in  relation  to  the  forward  Euler 
method  is  apparent  as  is  the  possibility  of  using  more  com- 
plex predictor-corrector  pairs  [62], 

Implicit  Integration 

Consider  now  an  alternative  approach.  The  approxima- 
tion to  the  derivative  x(r)  could  just  as  easily  have  been 
made  at  t + h.  In  this  case  (41)  becomes 

x(l  + A)  « A/(x(r  + A))  + x(t).  (45) 

This  is  known  as  the  backward  Euler  integration  formula. 
Since  x(r)  is  known,  (45)  represents  an  implicit  equation  in 
the  unknown  xfr  + A).  This  equation  may  be  solved  for 
x(t  + A)  by  the  Newton-Raphson  method  previously  con- 
sidered and  the  process  repeated  at  each  point  in  time.  This 
procedure  is  called  implicit  numerical  integration.  The 
trapezoidal  integration  formula  (42)  is  also  an  implicit 
integration  formula. 

All  of  the  methods  discussed  thus  far  are  of  the  general 
form 

k I 

t(t  + A)  » £ a,x(t  - iA)  + A £ bjXit  - jh).  (46) 
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It  is  convenient  to  introduce  an  alternate  notation.  Without 
loss  of  generality  it  can  be  assumed  that  the  time  step  ft  is 
constant  and  thus  that  l»nli.  The  notations!  change  is 
defined  such  that 

x(t  + ft)  « i((n  + lVi)  = jt.+ (47) 
Equation  (46)  can  now  be  written  as  (48) 

*.♦!  - I ‘V'.-i  + h X bji.-j.  (48) 

i-0  /*-! 

Methods  used  in  which  ft. , =0  such  as  (41)  or  (43)  and  (44) 
are  termed  explicit  integration  methods  while  methods  used 
in  which  ft.,  #0  such  as  (42)  or  (45)  are  termed  implicit 
integration  methods. 

Runqt-Kutla  Integration 

A third  class2  of  integration  methods  differs  from  the 
preceding  two  classes  in  that  the  interval  between  t and  f + ft 
is  divided  into  subintervals.  These  routines,  which  are 
referred  to  as  Runge-Kutta  methods  [43],  seek  to  establish  a 
very  good  approximation  to  an  average  derivative  in  the 
subinterval.  A fourth-order  Runge-Kutta  procedure  is 
given  by 


- 

Ml  + ft)  - 40  + g(/,  + 2 f2  + 2/, 

+ /a)  (49) 

where 

- 

/,  - f(xUlt) 

(50a) 

I,  -/(*!  * *5) 

(50b) 

ft  - /(mi)  + ^ 

(50c) 

/«  * /^*<0  + fr/i.  ‘ + 2) ' 

(50d) 

In  addition,  both  higher  and  lower  order 

Runge-Kutta 

methods  are  available. 

Stability  of  Numerical  Integration 
With  the  preceding  material  as  background,  the  problem 
of  step-size  determination  with  regard  to  stability  can  be 
considered  (IS).  In  this  context,  stability  refers  not  to  the 
response  of  the  physical  circuit  but  rather  to  whether  or  not 
the  errors  generated  at  each  step  of  the  numerical  integra- 
tion process  tend  to  decay  (stable)  or  grow  (unstable)  as 
the  solution  progresses  in  time. 

The  study  of  stability  is  carried  out  by  substituting  the 
differential  equation  of  interest  into  the  particular  form  of 
the  integration  formula  being  considered  and  generating  a 

1 An  alternative  classification  scheme  is  baaed  on  the  number  of  pre- 
vious time  points  at  which  values  ot‘  air)  or  air)  are  required  to  compute 
alr-r-A).  Single-step  methods  require  values  only  at  limes  greater  than  or 
* equal  tor.  Mullmep  methods  may  requ’ re  values  at  r -A  r-  2A.  etc.  In  this 

tense,  the  Runge-Kutta  method  described  here  is  a single-step  method. 
It  is  also  possible  to  use  Runge-Kutta  methods  in  a multutep  manner  by 
considering  subtntervals  of  width  h in  a total  interval  of  width  eh  where 
n>  I. 


difference  equation.  The  roots  of  the  difference  equation  are 
found  by  converting  the  difference  equation  into  an  alge- 
braic equation,  i In  some  cases  the  ; transform  can  be  used.) 
If  the  order  of  the  difference  equation  exceeds  the  order  of 
the  original  differential  equation,  here  assumed  to  be  first 
order,  then  all  but  one  of  the  solutions  to  the  difference 
equation  are  parasitic  Stability  requires  that  all  parasitic 
solutions  have  magnitudes  less  than  unity.  Those  parasitic 
solutions  whose  magnitudes  exceed  unity  represent  growing 
solution  modes  whtcs  can  be  excited  by  local  truncation 
errors  and  thus  dominate  the  correct  solution. 

Through  a transformation,  the  solutions  to  the  differ- 
ence equation  can  be  related  to  the  eigenvalues  of  the  circuit 
at  each  time  point  and  the  time  step  ft.  For  explicit  integra- 
tion and  Runge-Kutta  methods.1  the  time  step  ft  must 
inevitably  be  held  smaller  than  a constant  divided  by  the 
largest  eigenvalue  at  the  present  time  point  in  order  to  keep 
the  parasitic  solutions  small.  The  constant  involved  is  usu- 
ally less  than  ten  (43  ].  On  the  other  hand,  for  implicit  meth- 
ods. it  is  found  that  this  requirement  may  be  relaxed  by 
three  or  more  orders  of  magnitude  [63].  As  brought  out 
earlier,  an  iterative  method  such  as  Newton-Raphson  must 
be  used  to  determine  the  state  of  the  circuit  at  t+h.  The 
additional  computation  time  required  to  solve  the  nonlinear 
system  of  equations  at  a time  point  is  more  than  offset  by 
the  comparatively  large  steps  in  simulation  time  which  may 
be  taken  between  points.  Note  that  in  this  regard  both  the 
backward  Euler  and  the  trapezoidal  methods  are  found  to 
be  stable  for  all  positive  values  of  ft  when  the  eigenvalues 
lie  in  the  open  left  half-plane.  This  does  not  mean  that  cir- 
cuits such  as  oscillators  and  multivibrators  whose  eigenval- 
ues may  lie  in  the  right  half-plane  cannot  be  analyzed  but 
rather  that  the  maximum  value  which  ft  may  be  allowed  to 
assume  must  be  reduced. 

Truncation  Error 

Stability  of  a numerical  integration  method  implies  only 
that  parasitic  solutions  when  excited  will  not  grow  with 
time.  Thus  stability  only  guarantees  that  as  time  is  allowed 
to  go  infinity,  the  computed  solution  will  converge  to  the 
exact  solution.  At  finite  times  the  solutions  may  differ  sig- 
nificantly due  to  truncation  error.  The  error  term  asso- 
ciated with  the  forward  Euler  method  mentioned  previously 
is  an  example.  In  that  particular  case,  the  error  is  associated 
with  h2  multiplied  by  the  second  derivative  or  curvature  of 
x(t).  Thus  in  regions  where  the  response  is  rapidly  varying, 
ft  may  have  to  be  chosen  many  times  smaller  than  stability 
considerations  require.  The  time  step  ft  can  be  increased 
when  the  response  is  slowly  changing. 

In  general,  when  a higher  order  integration  formula  is 
used,  the  truncation  error  may  be  made  proportional  to 
higher  derivatives.  Similar  to  a Taylor  senes,  higher  order 
terms  tend  to  zero.  The  order  of  an  integration  formula  can 
be  determined  in  the  following  way.  If  linear  differential 
equations  whose  solutions  are  polynomials  of  finite  order 
are  considered,  the  order  of  the  polynomial  of  largest  degree 

1 Here  only  second  or  higher  order  methods  are  implied  where  order 
is  denned  in  the  subsection  on  truncatusn  error.  For  first-order  methods, 
no  parasitic  solutions  exist. 
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for  which  (48)  is  exact  is  the  order  of  the  method.  Without 
further  justification,  it  is  slated  that  backward  Euler  is  first 
order  while  trapezoidal  integration  is  second  order. 

A difficulty  with  higher  order  methods  is  that  the  number 
of  parasitic  solutions  to  the  stability  difference  equation  is 
increased  with  the  order,  and  the  restrictions  on  h become 
increasingly  severe.  A compromise  must  be  made  between 
using  highly  stable  methods  with  larger  truncation  errors 
and  less  stable  methods  with  reduced  truncation  errors.  An 
approach  used  by  several  writers  [64f-[66j  is  to  vary  the 
order  of  the  method  as  the  calculation  proceeds,  based  on 
an  examination  of  the  eigenvalues  or.  equivalently,  the  rate 
at  which  the  response  is  varying.  The  aim,  of  course,  is  to 
use  as  high  an  order  method  as  possible  consistent  with 
stability  [64], 

\onlinear  Transient  Analysis  Programs 

With  the  exception  of  the  Runge-ICutta  approach,  the 
numerical  integration  formulas  considered  thus  far  and 
generalized  in  terms  of  (48)  have  all  been  linear  in  the  sense 
that  the  value  xm+ , is  expressed  as  a linear  combination  of 
function  values  and  derivatives.  Pope  [67]  and  Fowler  and 
Warten  [68 1 have  developed  modifications  of  predictor- 
corrector  methods  in  which  an  exponential  term  is  included. 
This  approach  allows  a larger  step  size  to  be  used  in  com- 
parison with  straight  predictor-corrector  or  Runge-ICutta 
methods  and  has  been  used  in  both  sceptre  and  circus. 
The  approach  does  not  allow  the  dramatic  increase  in  time 
steps  which  is  possible  with  an  implicit  method,  sceptre 
includes  two  additional  integration  routines,  a trapezoidal 
predictor-corrector  method  similar  to  (43)  and  (44)  and  a 
fourth  order  Runge-ICutta  method  similar  to  (49)  and  (50). 

cutPAC.  though  not  generally  available  outside  of  Bell 
Telephone  Laboratories,  demonstrates  the  advantages  of 
implicit  integration  methods.  The  program  uses  the  second- 
order  implicit  integration  formula  [43] 

*.♦!*-  iJC—i  + T*.  + S64,a.|  (51) 

which  is  given  here  without  justification.  The  step  size  h, 
which  is  allowed  to  vary,  is  Wept  as  large  as  possible  consis- 
tent with  maintaining  a small  local  truncation  error. 
Shichman  [43)  reports  that  cirpac  typically  runs  up  to 
ten  times  faster  than  circus  and  up  to  twenty  times  faster 
than  an  earlier  version  of  cirpac  which  used  a predictor- 
corrector  routine. 

The  programs  considered  to  this  point  use  state-variable 
formulations  where  the  presence  of  a first-order  differential 
equation  is  made  readily  apparent.  The  trac  program 
which  employs  nodal  analysis  U3cs  a trapezoidal  implicit 
integration  method.  The  formulation  can  be  conveniently 
illustrated  for  the  case  of  a linear  capacitor  where  the  i-o 
characteristic  is  given  by 


In  integral  form.  (52)  can  be  rewritten  as 

J*r»»  (•>»» 

it(l)df  = C|  ii[u,(t)  “ 0j(()] 
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where  tre(f|  is  taken  to  be  the  difference  between  the  node 
voltages  r.Ul  and  Mr).  The  trapezoidal  integration  formula 
(42)  applied  to  the  lelt-liaud  side  ol  |53|  yields 

q [if(r  + (i)Tie(tl]  = C[r,(t  vhl— e,(r|  — CjU  + hl+DjUl].  154) 

This  equation  can  be  rewritten  to  ootauf 
’C 

ic(t  + h)  = ^-[PjU  +•  hi  - CjU  t hi] 

''C  . 

- [mi I - 11,(1 )]  - i,(r).  (55) 

This  Iasi  expression  can  be  represented  in  terms  of  an 
equivalent  conductance  2C,h  in  parallel  with  a current 
source  between  nodes  i and  j of  value  — (2C/li)[t),-(t)  — r,(r|] 
- ie|r).  Both  elements  are  readily  handled  in  terms  of  nodal 
analysis. 

Other  reactive  circuit  elements  are  handled  by  trac  in  a 
similar  manner.  The  modified  Newton- Raphson  method 
used  for  dc  analysis  and  previously  described  is  used  at  each 
time  point.  While  nearly  as  fast  as  cirpac.  trac  does  not 
incorporate  the  same  sort  of  variable  step-size  feature  and 
truncation  error  control.  Rather  the  user  has  the  option  of 
specifying  up  to  ten  time  intervals  and  the  maximum  step 
size  to  be  allowed  in  each  time  interval.  Thus,  where  it  is 
known  that  a switching  transient  is  about  to  occur,  the  user 
can  force  a smaller  step  size  to  be  imposed.  An  advantage 
of  trac  is  that  it  is  readily  available.  Further,  it  is  a relatively 
small  program  which  consists  of  only  2000  statements. 

Recently,  a modified  version  of  trac  known  as  tracap 
and  released  as  a part  of  the  syscap  package  has  been 
announced.  Like  dicap  its  dc  counterpart,  tracap  is  only 
available  to  users  of  Control  Data  Corporation  s remote 
batch  service  Cybernet.  A version  of  trac  known  as 
mtrac  which  handles  magnetic  cores  has  also  been  devel- 
oped [27], 

The  transient  portion  of  ecap  is  restricted  to  circuits  in 
which  nonlinear  elements  are  described  by  fixed  piecewise- 
dnear  models.  Ideal  switches  are  used  to  move  from  one 
linear  segment  of  a model  to  another  in  accordance  with 
the  direction  of  current  through  a sensing  branch.  The  use 
of  this  portion  of  the  program  is  quite  cumbersome  if 
piecewise-linear  reactive  elements  are  included.  Implicit 
integration  is  used. 

As  to  the  future,  net-II4  a completely  revised  version  of 
net-1  (in  fortran)  is  under  extensive  test. 

VI.  Linear  Pole-Zero  Analysis 

Linear  pole-zero  circuit  analysis  is  considered  as  a sep- 
arate topic  because  of  the  number  of  different  techniques 
used  and  the  specialized  problems  involved  relative  to  fre- 
quency-domain analysis.  Naturally,  the  poles  and  zeros  of 
the  transfer  function,  once  obtained,  can  provide  the  same 
frequency  response  information  as  (he  linear  ac  analysis 
programs  previously  considered.  However,  in  many  prob- 
lems it  is  the  poles  and  zeros  themselves  which  are  of  prime 
interest  to  the  circuit  designer. 


(53) 


VET-II  is  J new  version  oinet-I  by  A.  Malm  Derg,  now  being  developed. 
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State- Variable  and  Eigenvalue  Approach 

The  cornap  program  wntten  by  Pottle  has  found  wide- 
spread use.  As  mentioned  previously  it  employs  a state- 
variable  formulation  leading  to  a set  of  equations  of  the 
form  of  (3)  and  (4).  In  the  complex  frequency  domain  these 
equations  take  the  form 

sx  = Ax  + Bu  (56a) 

y = Cx  + Du.  (56b) 

The  output  vector  y becomes 

y = [C(sf  — A)~lB  + D]«.  (57) 

For  the  single-input-single-output  case,  the  circuit  transfer 
function  T(s)  is  given  by 

T(s)  = = Cisl  - Ar'B  -i-  D (58) 

u(s) 

where  D is  a scalar.  The  poles  of  the  transfer  function 
(natural  frequencies  of  the  circuit)  are  given  by  the  zeros  of 
det  (si— A). 

In  an  early  version  of  cornap  the  recursive  method  of 
Levemer  [69],  also  known  as  the  Souriau-Frame  algo- 
rithm [70],  was  used  to  construct  the  adjoint  of  (si -A) 
and  the  characteristic  polynomial  associated  with  A. 
Muller’s  method  (71  ] was  then  used  to  find  the  zeros  of  the 
polynomial.  The  recursive  method  however  was  particu- 
larly sensitive  to  roundoff  errors  which  severely  limited 
accuracy. 

In  the  present  version,  use  is  made  of  the  fact  that  the 
zeros  of  det  (si  - A ),  which  are  the  zeros  of  the  characteristic 
polynomial  and  the  poles  of  T(s).  are  also  the  eigenvalues 
of  the  matrix  A.  The  QR  algorithm  of  Francis  [72J-[74]  is 
used  to  compute  these  eigenvalues. 

The  transmission  zeros  of  a circuit  transfer  function  are 
obtained  by  making  use  of  the  fact  that  the  zeros  of  the 
transmission  function  of  a feedback  network  placed  around 
an  ideal  amplifier  of  infinite  gain  are  the  poles  of  the  closed- 
loop  transfer  function.  The  transfer  function  T(s)  is  con- 
sidered to  be  placed  around  such  an  ideal  amplifier.  State 
equations  describing  this  inverse  system  [75]  are  obtained 
in  terms  of  the  original  A,  B.  C.  and  D matrices.  The  eigen- 
values of  the  inverse  system  are  then  the  zeros  of  the  original 
system. 

cornap  can  handle  up  to  32  state  variables  (eigenvalues) 
but  does  require  a 64  x 64  double  precision  matrix  plus  an 
additional  32  x 64  double  precision  matrix.  The  core  re- 
quirements are  thus  large  and  it  has  been  noted  that  accu- 
racy is  questionable  for  some  large  circuits  of  the  order  of 
64  branches  and  24  nodes. 

A recent  paper  by  Sandberg  and  So  [76]  describes  an 
alternate  approach  for  obtaining  transmission  zeros.  The 
approach  still  makes  use  of  the  QR  algorithm  for  finding 
eigenvalues. 

Iterative  Technique 

The  trfn  portion  of  lisa  uses  a very  accurate  iterative 
procedure  to  find  the  natural  frequencies  of  a circuit.  A 
nodal  formulation  is  used  and  each  term  of  the  nodal 


admittance  matrix  represents  a polynomial  in  s.  the  complex 
frequency  variable.  For  single  input  single  output,  the 
transfer  function  is 


, , r.(s)  Ori 

7. I s ) — s ■ 

i.(s)  det  Ms) 


(59) 


T„.(s)  is  the  minor  of  the  element  y„(s),  the  poles  of 
7[«.(j)  are  the  zeros  of  det  T(j).  The  zeros  of  T_u)  are  the 
zeros  of  det  ^(j).  The  iterative  root-finding  method  of 
Muller  is  applied  directly  to  det  f(.vi  and  det  l(«.(j|  In 
Muller’s  method,  the  determinant  is  initially  evaluated  at 
three  points  and  modeled  by  a quadratic.  A zero  of  this 
quadratic  is  then  used  in  place  of  one  of  the  original  evalua- 
tion points.  This  procedure  is  continued  and  a sequence  of 
better  and  better  approximations  to  each  zero  of  the  deter- 
minant is  developed.  The  amount  of  determinantal  evalua- 
tion time  appears  excessive;  however,  the  accuracy  is 
excellent.  Further,  it  is  possible  to  program  this  approach 
very  effectively  to  reduce  greatly  the  computation  time. 
frank,5  a program  of  this  type,  is  faster  for  large  circuits 
than  cornap,  more  accurate  requiring  only  single  precision, 
and  uses  significantly  less  core.  Finally,  where  element 
values  are  to  be  varied  and  analyses  repeated,  a large  savings 
in  time  may  result  because  previous  solutions  can  be  used 
to  supply  good  initial  estimates  of  new  solutions. 


Laplace  Expansion  Approach 

By  a straightforward  Laplace  expansion  of  the  deter- 
minant of  the  nodal  admittance  matrix  T(s).  the  coeffi- 
cients of  the  characteristic  polynomial  can  be  obtained. 
Muller's  method  can  then  be  used  to  solve  for  the  zeros  of 
the  polynomial.  An  approach  of  this  type  is  used  in  the 
poly  portion  of  lisa;  however,  it  is  severely  limited  by  the 
effects  of  roundoff  error  and  is  limited  to  circuits  of  less 
than  12  nodes. 

Larger  circuits  can  be  handled  effectively  with  programs 
of  this  type  if  a restriction  is  made  to  active  RC  circuits. 
Several  polynomial  manipulations  are  then  eliminated  and 
accuracy  and  speed  are  improved.  Program  sprague"  has 
been  found  to  be  accurate  and  fast  for  circuits  up  to  19 
nodes. 


Nodal  Analysis-Eigenvalue  Approach 
The  eigenvalue  approach  based  on  a nodal  equation 
formulation  can  be  used  for  circuits  restricted  to  active  RC 
elements.  The  nodal  admittance  matrix  has  the  form 

JTs)  - G + sC  (60) 

Eigenvalue  programs  of  this  type  have  been  developed  both 
at  M.I.T.  [30],  [31]  and  at  the  Technical  University  of 
Denmark  [77],  The  M.I.T.  program  requires  the  user  to 
enter  the  elements  of  the  matrices  G and  C of  (60)  directly  ; 
however,  a topological  input  routine  can  easily  be  added. 

’ frank  is  a program  written  and  developed  in  1969  6y  the  1C  Group. 
Electronics  Research  Laboratory. -University  of  California.  Berkeley 
• Sprague  is  a program  wntten  and  developed  in  I96S  by  the  1C 
Group.  Electronics  Research  Laboratory.  University  of  California. 
Berkeley. 
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The  method  is  based  on  the  transformation  of  det 
(C-f  jC)  to  del  i(jC~  ^ — sl)  where  / is  the  identity  matrix 
and  6 and  C~  ' are  matrices  of  rank  less  than  or  equal  to  n 
the  order  of  Vis).  The  transformation  is  earned  out  such 
that  the  zeros  of  the  original  determinant  are  unchanged. 
The  zeros  of  the  alternate  determinant,  however,  are  also  the 
eigenvalues  of  6 C~'.  The  QR  method  is  used  to  compute 
these  eigenvalues.  The  transmission  zeros  are  obtained  by 
applying  the  same  procedure  to  n-isi- 
If  the  rank  of  C is  equal  to  i,  6C~ 1 can  be  obtained  by 
applying  a variation  of  Gaussian  elimination  known  as 
Gauss-Jordan  reduction  to  C and  performing  the  same 
operation  on  G.  Gauss-Jordan  reduction  differs  from 
Gaussian  elimination  in  that  a matrix  is  diagonalized  and 
then  normalized  to  obtain  the  identity  matrix.  This  requires 
elimination  of  elements  above  the  diagonal  as  well  as  those 
below.  If  the  rank  of  C is  less  than  n.  then  redundant  rows 
and  columns  exist  in  C and  C and  must  be  eliminated  (77  ]. 

This  method  again  requires  only  sufficient  core  to  store 
the  nodal  admittance  matrix  and  is  very  last.  However,  it  is 
a relatively  new  method  and  its  relative  accuracy  has  yet 
to  be  determined,  though  it  should  be  better  than  that  of  the 
Laplace  expansion  approach. 

An  extension  to  active  RLC  networks  has  also  been  de- 
scribed [78],  Inductive  elements  are  replaced  by  ideal 
gyrator-capacitor  equivalent  circuits.  The  effectiveness  of 
this  extension  to  the  preceding  method  has  not  yet  been 
evaluated,  however. 

Summary 

A number  of  techniques  which  are  used  for  computing 
the  poles  and  zeros  of  circuit  transfer  functions  have  been 
considered.  Though  each  approach  requires  more  extensive 
computations  than  are  required  to  compute  responses  in 
the  frequency  domain,  the  additional  insight  obtained  from 
a knowledge  of  the  poles  and  zeros  often  merits  the  addi- 
tional effort.  In  all  cases,  the  order  of  the  circuits  which  can 
be  analyzed  is  significantly  less  than  that  which  can  be 
handled  with  the  frequency-domain  program  described 
earlier. 

No  mention  has  been  made  of  programs  using  a topolog- 
ical tree  enumeration  formulation.  While  accurate  (79  J- 
(80 1 and  useful  for  sensitivity  studies  because  of  the  explicit 
form  of  the  coefficients  of  the  polynomials,  these  programs 
have  been  found  to  be  too  slow  to  be  of  interest  and  are 
limited  to  very  small  circuits.  Similar  conclusions  have  been 
found  to  apply  to  programs  based  on  how  graph  methods 
[81). 

VII.  Conclusion 

In  the  preceding  sections  of  this  paper,  the  basic  elements 
of  computer-aided  circuit  analysis  have  been  reviewed  with 
emphasis  on  those  techniques  and  routines  necessary  for 
the  adequate  simulation  of  four  basic  classes  of  circuits: 
linear  dc  and  ac.  nonlinear  dc.  nonlinear  transient,  and 
linear  pole  zero.  One  topic  not  considered  but  very  impor- 
tant is  the  development  of  device  models  suitable  for 
computer-aided  circuit  analysis.  An  overview  of  modeling 
would,  * itself,  require  a full  paper  and  thus  cannot  be 


:5 


undertaken  here.  However,  several  articles  on  Ihe  modeling 
of  bipolar  and  lidd-eflect  transistor  models  are  included 
in  (82|-185|. 
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Equations  for  a Sparse  Matrix  Solution 
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Abstract — Spars a matrix  atorapa  and  tohition  tacbmquaa  ara  used 
iitaMivdy  in  vary  Iqrgq  ayatama  of  hitndrada  of  linear  aqua- 

tions wfwcb  arias  in  the  analyses  of  muitipiy  intarconnsetsd  pfiyaecai 
■yatama.  Thaaa  taebniquaa  hays  oft  an  baan  ovarloobad  in  the  analysts 
of  raiatnrofy  small  aiactric  networks  aeon  though  thoir  use  can  result 
in  vary  significant  improvements  in  computer  storage  requirements 
and  sxecutton  times.  The  time  savings  is  particularly  noticeable  whan 
many  solutions  for  the  same  circuit  with  different  parameter  values 
ara  required. 

A particular  sparse  matria  storage,  reordering,  and  solution  tech- 
nique is  d ascribed.  A node  renumbering  algorithm  which  is  specifi- 
cally directed  at  preserving  the  sparse  structure  of  nodal  admittance 
matrices  during  the  solution  by  Gaussian  elimination  is  described  in 
detail.  Computer  flow  charts  for  the  renumbering  ara  included  along 
with  specific  circuit  seam  plea  which  compare  the  relative  computa- 
tional effort  required  for  spares  solution  versus  full  matrix  solution. 

I.  Introduction 

MANY  circuit  analysis  programs  that  are  in  wide- 
spread use  today  are  limited  by  the  necessary 
' storage  requirements  of  the  admittance  or  im- 
pedance matrix  and  are  inefficient  in  obtaining  the  solution 
to  the  set  of  equations  describing  the  network.  If  nodal 
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analysis  is  used,  the  resulting  admittance  matrix  of  the 
circuit  is  virtually  always  sparse ; that  is,  less  than  50  percent 
of  the  matrix  elements  are  nonzero. 

Many  authors  have  written  on  the  subject  of  solutions  of 
sparse  systems  of  equations  [I  ]—  [3  ].  Results  have  indicated 
that  time  savings  of  orders  of  magnitude  can  be  achieved 
for  solution  of  large  networks.  This  paper  describes  a solu- 
tion algorithm  which  was  specifically  developed  for  nodal 
analysis  of  electric  circuits.  The  algorithm  takes  advantage 
of  the  sparse  nature  of  typical  network  admittance  matrices 
to  obtain  the  voltage  solution  vector  by  Gaussian  elimina- 
tion. Storage  is  allocated  for  only  nonzero  terms  of  the 
admittance  matrix  and  Gaussian  elimination  proceeds 
through  only  these  terms,  thereby  allowing  efficient  analysis 
of  much  larger  circuits  than  was  previously  possible. 

A fortran  iv  program  utilizing  the  sparse  storage  and 
solution  algorithms  has  been  written  and  is  currently  in  use 
in  a linear  ac  analysis  program.  It  can  handle  100  node 
circuits  in  less  than  30  000  words  of  storage  on  the  campus 
CDC  6400  computer  1 The  29-node  circuit  shown  in  Appen- 

' Thu  program  was  developed  by  graduate  students  in  the  Department 
of  Engineering.  University  of  California.  Berkeley  during  1909-1970, 
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dix  II  required  just  3.5  s to  obtain  71  complete  solutions  at  a 
like  number  of  frequency  points.  The  time  includes  the  re- 
loading of  the  admittance  matrix  for  each  of  the  71  fre- 
quency points. 

The  paper  is  divided  into  three  sections  and  has  Appen- 
dixes. Section  II  describes  the  system  of  pointers  developed 
for  rapid  access  of  the  nonzero  elements  of  any  network. 
Section  III  describes  the  decomposition  and  optimum  re- 
numbering scheme  to  preserve  the  sparse  matrix  structure 
during  Gaussian  elimination  [4 1 on  the  admittance  matrix. 
Section  IV  describes  the  decomposition  and  solution  algo- 
rithms which  operate  only  upon  the  nonzero  terms  of  the 
matrix.  Appendix  I includes  the  flow  charts  for  the  renum- 
bering scheme.  Appendix  II  illustrates  two  circuit  examples. 

II.  Admittance  Matrix  Storage 

The  admittance  matrix  Y for  a typical  ten-node  electric 
circuit  contains  fewer  than  50-percent  nonzero  terms.  As  the 
number  of  circuit  nodes  increases  the  percentage  of  non- 
zero terms  drops,  so  that  at  100  nodes  typical  circuits  con- 
tain 5-percent  nonzero  terms.  For  efficient  utilization  of 
high-speed  memory  and  to  allow  for  practical  solution  of 
very  large  circuits,  storage  is  allocated  for  only  the  nonzero 
terms  of  the  admittance  matrix.  These  terms  are  collapsed 
into  the  columnar  array  A.  An  efficient  set  of  pointers  for 
locating  these  terms  in  the  array  is  an  absolute  necessity. 

The  matrix  is  symmetric  except  for  terms  associated  with 
solid-state  device  equivalent  circuits.  For  the  most  part, 
these  terms  are  structurally  symmetric  but  numerically  non- 
symmetric.  Pointers  locate  nonzero  terms  and  therefore  can 
take  advantage  of  structural  symmetry,  even  though  nu- 
meric symmetry  does  not  exist. 

The  columnar  array  A of  nonzero  terms  is  organized  into 
the  three  sections  shown  in  the  following  table. 


The  first  section  ofthe  4 array  is  reserved  All) «>,, 

lor  the  diaeonai  term*  of  the  admittance  4(2)-y„ 
auim  1 1 Mellowed).  Position  * contains  4(3)-v„ 
the  vM  diagonal  term;  [hereto re,  no 
portlier  rynwn  u necessary 


The  pointer  system  for  these  last  two  sections  takes  ad- 
vantage of  the  symmetric  structure,  so  that  one  set  of  point- 
ers serves  both  of  the  triangular  portions. 

The  pointers  are  set  up  in  two  stages.  In  the  first  stage,  the 
nonzero  structure  of  the  nodal  admittance  matrix  is  re- 
corded by  the  pointer  system  as  the  individual  network  ele- 
ments and  their  respective  node  connections  are  transferred 
from  data  cards  to  the  program.  Based  upon  this  structure, 
a node  renumbering  is  effected  which  attempts  to  minimize 
the  number  of  operations  required  for  triangular  decompo- 
sition of  Y by  Gaussian  elimination.  After  completing  the 
renumbering,  the  pointer  system  for  the  lower  and  upper 
triangular  portion  of  Y is  changed  to  correspond  to  the  new 
node  numbers.  Storage  of  admittance  terms  in  the  second 
and  third  sections  of  the  columnar  array  A is  to  be  based 
upon  the  renumbered  system,  while  the  diagonal  section  is 
to  be  arranged  according  to  the  original  node  numbers. 
Only  the  pointers  associated  with  the  upper  triangular  array 
of  the  renumbered  system  of  equations  are  retained.  Num- 
bers of  the  lower  triangular  portion  can  be  located  by  using 
symmetry  and  the  upper  pointers. 

The  pointers  consist  of  two  integer  arrays.  The  first  array 
iur  contains  ,V  integer  numbers  where  .V  is  the  number  of 
circuit  nodes.  The  number  stored  in  position  k of  this  array 
represents  the  starting  location  in  the  second  pointer  array 
iuc  of  terms  associated  with  row  k of  the  admittance  matrix. 
In  stage  one  the  second  array  includes  the  column  location 
of  all  nonzero  olf-diagonal  terms  in  each  row  of  the  original 
Y matrix.  In  stage  two  the  positions  in  the  first  pointer  array 
and  all  of  the  numbers  of  the  second  pointer  array  corre- 
spond to  the  new  node  numbers.  Every  nonzero  term  that 
will  occur  in  the  final  decomposed  form  of  the  upper  tri- 
angular matrix  is  identified  by  this  final  pointer  set.  At  this 
stage  the  position  in  the  second  pointer  array  directly  corre- 
sponds to  the  position  in  the  collapsed  admittance  array  of 
the  term  identified.  As  an  example  to  illustrate  the  pointer 
system,  consider  an  admittance  matrix  with  the  following 
pattern  of  nonzero  terms  before  node  renumbering,  where 
the  y, j are  the  nonzero  terms: 


The  second  section  of  the  A array  is  re- 
served for  the  nonzero  off-diagonal  terms 
of  the  tipper  tnenguter  portion  of  the  ad- 
mittance matrix  (400  totRl  allowed). 
These  terms  are  stored  by  rows. 


A(it)mym 


>ll 

0 

y u 

0 

yis' 

0 

Yil 

yu 

yw 

0 

4<IM) 

Yu 

•I’ll 

yjs 

yia 

0 

4(10l)-y,J 

0 

ya  2 

yal 

yea 

0 

4(l02)-y„ 

.Vst 

0 

0 

0 

y5s. 

The  pointer  arrays  before  node  renumbering  would  con- 
tain the  following  numbers. 


The  third  section  of  the  A array  is  re-  4(500) 
served  for  the  nonzero  off-diagonal  terras  4(50l)-vM 
of  the  lower  triangular  ponton  of  the  ad-  4(102)- v,, 
mitlancc  matrix  (400  total  allowed). 

These  terms  are  stored  by  columns. 

4(/n  + 400)-.y,..l 


4(900) 


Row  Locator 

Column  identifier 

Term  Identified 

ium<n-  l 

tuct  1 > »3 

Y\s 

run(2)-3 

iua2)  -5 

fts 

iu*(3)»5 

iua3)  »3 

IUC14)  »4 

y:* 

10 

iuc<5)  »l 

y»» 

iuw(6)-ll 

kjc(6)  »2 

r» 

JUCI7)  »4 

iuc(8)  -2 

y 4s 

tuc(9l  -3 

V.  J 

iuc(10)*  | 

y si 
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1 

on  the  right. 

the  following: 

L 

u 

m 

r 

1 

0 

0 

O' 

“11 

“12  “l 3 

“14' 

■>u 

>12  >11 

>14 

'2, 

1 

0 

0 

0 

“22  “22 

“24 

>21 

>22  >23 

>14 

'3. 

*32 

1 

0 

0 

0 “33 

“14 

>31 

>32  >11 

>14 

<4. 

*42 

*41 

1. 

0 

0 0 

“44 

>41 

>42  >41 

>44 

After  node  renumbering  the  admittance  matrix  might 
assume  the  following  form: 


The  matrix  U is  an  upper  triangular  matrix: 


>11 

0 

0 

>14 

>11 

0 

>12 

>23 

>24 

0 

0 

>32 

>31 

>34 

0 

>41 

>42 

>43 

>44 

>43 

.>11 

0 

0 

>34 

>33 

“12 

“22 

0 

0 


“13  “1* 

“23  “24 

“33  “34 

0 U«4 


The  numbers  within  the  pointer  arrays  would  be  the  fol- 
lowing with  the  admittance  terms  in  the  A array  as  shown 


A full  matrix  decomposition  will  be  illustrated  to  show  the 
steps  involved  in  a triangular  decomposition  by  Gaussian 
elimination.  The  product  of  the  L and  U matrices  shown  is 


IUR(I)-I 

iuctl)-a 

4(IOI)-y„ 

iu«(2)-3 

tuc(2)-5 

4(102)— yu 

iur(3)-5 

iUC(3)-3 

4(l03)-y,j 

iur(4)»6 

iuc(4)-4 

4(104t-yu 

nm<5)»7 

IUC(5)-4 

4(105)— y)4 

WC(6)-5 

4<l06>-y., 

The  Am  pouter  imy  ojr(*) 

The  seoood  pouter  array 

4(50l)-y., 

indicates  the  starting  loca- 

nxtj)  indicates  the  column 

4(502)-y„ 

uom  in  iv&J)  of  the  nonzero 

location  of  the  nonzero  term 

4(503)-y„ 

terms  identified  in  row  (cot- 

uni  in  4 ( 100  * /l  sod  the 

4(504)-y41 

uma)  it  of  the  admittance 

row  location  of  the  nonzero 

4(505)— y4J 

ntecru. 

tern  stored  in  ^(500-*-;). 

4<S06)-y,4 

By  examination 


*u 

“12 


1 I'll 

■ yi2 

' >13 
■>14- 


By  multiplying  the  rows  of  L by  the  first  column  of  U,  we 
see  that 


111.  OecoMrosmoN  and  Node  Renumbering 
The  first  step  toward  obtaining  the  solution  vector  x in  the 


*2l“l 1 

= >21 

or 

*21 

= >2l/“l 1 

*3 1 “ 1 1 

” >31 

or 

*31 

* >3  l/“l 1 

*4l“ll 

« >41 

or 

*4. 

“ >4 1/“  1 1 

set  of  equations  Yxmb  involves  triangular  decomposition 
of  Y,  the  symmetrically  structured  admittance  matrix,  by 
Gaussian  elimination.  This  technique  involves  the  breaking 
down  of  Y into  a product  of  two  unique  matrices,  L and 
V,  where  Y m LU.  The  matrix  L is  lower  triangular  in  form 
with  1 as  the  value  of  every  diagonal  element. 


The  remaining  u,;  and  lt)  values  in  terms  of  yi(  and  pre- 
viously determined  lA,  uk/,  and  uu  values  are 


*2l“l2  + “22  " >22  OT  “22  =*  >22  ~ *2l“l2 


1 


0 0 
1 0 
<32  1 


(4,  *. 
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*2l“l3  + “23  ~ >23  0r  “23 
*2l“l4  + “24  = >24  Of  “24 

*Sl“l2  + *ll“22  “ >32  °r  *32 


*4l“l2  + <42“32  ” >42  Or  Ul  - 


>23  ~ *2l“l3 
>24  ~ *2l“l4 
>32  ~ *3l“l2 
“22 

>42 •-  *4l“l2 
“22 


*ll“ll  + *32u23  + “33  “ >33  0r  “l3  = >33  ~ *3l“l3  ~ *3l“23 
*St“l4  + *3l“24  + “14  » >14  °f  «34  = >14  “ *3l“t4  ~ *3l“24 

'.,“,3  + *4l“l3  + *.,“,3  - >43  Or  U,  » 

U33 

*4»“ie  + Ui*74  + **3«34  + * y44  Of  tt44  - Yaa  - /4|“t4  “ UlU2A 


UiuS 


... 
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The  first  row  of  U and  the  first  column  of  L are  obtained 
quite  easily.  In  programming  it  is  very  desirable  to  replace 
the  values  of  the  yi(  with  the  values  of  either  the  u„  or  the  ltl 
as  they  are  determined.  Also,  as  each  /,,  is  determined  the 
remaining  values  of  ya  in  the  ith  row  for  k > I are  changed 
by  subtracting  from  them  the  product  I, ,uit.  After  this 
first  row  and  first  column  step  of  the  decomposition,  the 
values  in  the  Y storage  locations  will  be  the  following: 


u,,  ut2  u13  u,» 

'a  I A'i  Ai  A'i 

<*.  A'i  A'i  A'i 

f*,  Aii  A'i  Aii 

where 

A'i  - yjj  - 'ii“u  - “22 

A'i  “ y2J  - <2l“l3  - “23 

yi"  - y2«  - '2i“ia  - “2« 

A'i  m y32  — *3  l“l  2 
A'i  - y33  - *3l“l3 
A"  - y34  - <3i“i« 

AH  * y«2  - (4t“i2 

yii  * y»3  - f»i“t3 

ylu  - y«  - 


Note  that  the  values  Ai>  Ai>  and  Ai  are  equal  to  u22, 
Ujj,  and  u24,  respectively.  /]2  and  lt2  are  obtained  from 
Ai  and  Ai,  respectively,  by  simply  dividing  by  the  pivot 
element  u22: 


*32 

'*2 


“22 

M 

“22 


Just  as  in  the  first  row-column  decomposition  step,  after 
each  ll2  is  obtained  the  remaining  values  of  u„  in  the  ith 
row  for  k >2  are  changed  by  subtracting  the  product  il2u2t 
from  them.  The  values  now  found  in  the  Y storage  array 
are  the  following : 


“n 
*21 
*31 

Ui 

where 

A'i  - A'i  ~ <32“23 
A'i  - A'i  - *32“24 
Ai  - Ai  - *42“23 
Aii  ~ Aii  ~ * 42“24- 


“l  2 

“13 

“14 

“22 

“23 

“24 

*32 

Ai 

Ai 

*4, 

Ai 

AH 

The  decomposition  proceeds  by  row-column  steps  until 
all  the  elements  of  V and  L are  determined  and  replace  the 
elements  of  Y in  storage.  In  general,  at  the  ith  row-column 
step,  the  ith  row  values  of  U are  determined  from  the  equa- 
tion 

= AT  for  all  j 2 i 

and  the  ith  column  values  of  L are  determined  from  the 
equation 

yi- n 

/ j,  =»  — — • for  all  j > i. 

“« 

This  technique  is  associated  with  Doolittle  [4]. 

The  sparse  matrix  decomposition  uses  the  Doolittle 
method  but  only  carries  out  the  operations  that  result  in  a 
finite  change  in  an  element  value  stored  in  the  Y array. 
A finite  change  in  a stored  value  occurs  under  just  two  cir- 
cumstances. The  first  arises  when  solving  for  the  1#  term 
just  mentioned  and  occurs  whenever  the  yj}- " term  is 
nonzero.  (The  diagonal  term  u:J  is  always  finite  in  the  ad- 
mittance matrix).  The  second  circumstance  arises  just  after 
the  determination  of  a finite  l,,  term  when  all  of  the  y£_  11 
terms  for  k > i are  changed  to  y'ji  by  the  relationship 

y£->V ■ 

Once  Ij,  is  determined  to  be  nonzero,  a change  between 
yJJ  and  yj f1*  will  occur  only  when  u12  is  also  nonzero.  A 
node  renumbering  scheme  is  effected  which  attempts  to 
organize  the  nonzero  terms  of  the  original  admittance 
matrix  in  such  a fashion  as  to  cause  the  second  circumstance 
to  arise  only  when  y%~ 11  is  already  nonzero.  Otherwise  an 
additional  nonzero  entry  will  occur  in  the  matrix  and  will 
increase  the  total  computational  effort  required  for  the 
decomposition. 

There  are  three  basic  parts  to  the  renumbering  algorithm. 
All  parts  search  the  nonzero  structure  recorded  by  the  stage 
one  pointer  system  described  previously.  The  algorithm 
takes  advantage  of  structural  symmetry  whenever  possible 
to  speed  up  the  renumbering  process.  An  array  numoft  is 
set  up  which  records  the  total  number  of  nonzero  olf- 
diagonal  terms  associated  with  each  node.  numoff( k ) 
equals  the  total  number  of  these  terms  that  would  appear  in 
the  Y array  in  row  k. 

Part  I of  the  algorithm  searches  this  array  once  to  see  if 
there  are  any  nodes  with  only  one  nonzero  off-diagonal 
term.  If  one  is  found,  it  is  numbered  I and  the  array  ni/moff 
is  altered.  The  olT-diagonal  term  of  this  new  node  1 will  be 
located  in  some  column  j.  Because  of  the  symmetric  struc- 
ture, there  will  also  be  an  off-diagonal  term  in  row  j, 
column  1.  During  the  course  of  the  decomposition,  these 
two  mentioned  off-diagonal  terms  will  become  u,t  and  l,,. 
respectively.  The  only  yji’  term  that  will  be  altered  during 
the  first  row-column  step  of  the  decomposition  will  be  the 
yIt  term,  which  will  change  to 

A'i  * yjJ  ~ Vi/- 


1 


J 
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All  of  the  other  utt  and  lkl  for  k > 1 terms  are  zero. 

Since  all  diagonal  terms,  i.e..  all  yu,  are  always  non- 
zero to  start  with,  there  will  be  no  additional  nonzero  terms 
created  at  this  step  of  the  decomposition.  In  searching  for 
a node  to  number  2.  we  can  remove  from  consideration 
both  of  the  nonzero  terms  ut/  and  Neither  of  these  terms 
will  be  used  again  during  the  course  of  the  decomposition. 
To  simulate  their  removal  the  array  numoff  is  changed  by 
reducing  the  recorded  number  of  off-diagonal  terms  associ- 
ated with  node  j by  1 . If.  by  this  reduction  of  l . the  effective 
number  of  otf-diagonal  terms  associated  with  node  i is  now 
one  or  fewer,  then  the  yth  node  will  be  numbered  2 and  the 
process  repeated.  A single  sweep  through  the  array  numoff 
will  rapidly  pick  off  every  node  that  has  only  one  or  fewer 
effective  off-diagonal  terms.  Decomposition  of  these  single 
x off-diagonal  term  equations  will  cause  no  new  nonzero 

terms  in  the  matrix. 

A typical  example  showing  the  working  array  after  the 
first  row-column  decomposition  step  is  shown  in  the  follow- 
ing. Row  one  has  just  one  nonzero  off-diagonal  term.  The 
array  numoff  contains  the  effective  number  of  nonzero  off- 
diagonal  terms  shown  on  the  right  after  1 is  subtracted  from 
row  5 : 


“ll 

0 

0 

0 

“is 

0 

NUMOFFI 1)  = l 

0 

y« 

0 

yj» 

0 

y« 

numoff<2)  - 2 

0 

0 

yjj 

0 

yss 

0 

NUMOFFI  3)  =»  1 

0 

y« 

0 

y*a 

y«s 

0 

NUMOFFI4)  =■  2 

is. 

0 

yss 

ysa 

yss 

0 

NUMOFFl  S)  » 2 

.0 

y«> 

0 

0 

0 

y»*. 

numoff(6)  =•  1. 

Part  II  of  the  algorithm  searches  the  remaining  nodes 
(those  not  renumbered  in  Part  I)  for  nodes  which  can  be 
decomposed  without  increasing  the  number  of  nonzero 
terms.  Suppose  m—  1 nodes  were  renumbered  by  the  Part  I 
search  and  suppose  node  i has  associated  with  it  two  non- 
zero off-diagonal  terms  (numoff(i)— 2)  in  the  portion  of 
the  matrix  that  has  not  been  renumbered  yet.  Let  these 
terms  be  ytl  and  y„.  Should  we  choose  to  renumber  this 
node  as  node  m,  the  decomposed  terms  would  become 
umt  and  u— . By  symmetry  there  would  also  be  a calculation 
of  the  finite  terms  /*,  and  I*,.  Two  terms  in  row  j and  two 
terms  in  row  k would  be  altered  by  the  decomposition. 
These  terms  are 

C - 

/T-ylT"-1*"- 
ytr-yfr1’ 
jC-jfi-'*-  u,. 

The  first  and  fourth  terms  are  on  the  diagonal  and  are  al- 
ready nonzero.  The  second  and  third  terms  are  off-diagonal 
but  are  in  symmetric  locations.  Therefore,  only  one  of  these 
terms  need  be  checked  to  see  if  it  is  nonzero.  If/J"  " should 
happen  to  be  nonzero,  then  the  node  would  be  renumbered 


next  and  no  change  in  the  structure  of  the  Y array  would 
occur  during  the  decomposition.  The  j and  k row  would 
have  1 removed  from  the  effective  number  of  off-diagonal 
terms,  and  if  this  caused  the  effective  number  to  drop  to  I in 
either  row  j or  k,  then  that  particular  row  would  be  re- 
numbered next. 

As  each  node  is  checked,  an  array  ifill  is  set  up  which  re- 
cords the  number  of  new  positions  that  would  become  non- 
zero if  that  particular  node  were  renumbered  next. 

After  checking  all  of  the  nonrenumbered  nodes  for  those 
which  can  be  decomposed  without  causing  an  increase  in 
the  nonzero  terms,  a check  is  made  to  see  if  any  were  re- 
numbered before  the  last  entry  into  the  Part  II  algorithm. 
If  any  nodes  were  renumbered,  the  algorithm  is  repeated 
because  now  the  effective  number  of  nonzero  off-diagonal 
terms  is  different  from  the  time  Part  11  was  first  entered.  It  is 
still  possible  that  additional  nodes  can  be  renumbered  that 
will  not  increase  the  nonzero  terms  during  decomposition. 
When  a complete  Part  11  search  is  made  without  finding  any 
nodes  tor  renumbering,  then  Part  111  is  entered.  At  this 
point,  every  remaining  node  would  cause  a change  in  the 
nonzero  structure  if  it  were  to  be  renumbered  next.  The 
array  ifill  indicates  how  many  new  nonzero  terms  would 
be  created  during  the  decomposition  if  the  node  in  question 
were  to  be  renumbered  next. 

Part  III  finds  the  node  that  would  cause  the  fewest  new 
nonzero  terms  by  searching  the  array  ifill.  In  case  more 
than  one  node  would  cause  the  same  fewest  new  nonzero 
terms,  the  node  among  these  with  the  most  number  of 
effective  nonzero  off-diagonal  terms  recorded  in  numoff  is 
numbered  next.  The  reason  for  this  choice  is  that  the  number 
of  new  nonzero  terms  created  during  the  decomposition 
at  this  step  is  still  a minimum,  but  also  the  number  of  terms 
removed  from  the  nonrenumbered  portion  is  the  largest 
possible,  subject  to  the  first  constraint.  By  this  choice  the 
sum  of  the  numbers  in  the  array  numoff  is  minimized. 

After  the  choice  is  made  and  a node  renumbered,  the 
new  nonzero  topology  caused  by  decomposition  of  that 
nodal  equation  is  recorded  in  the  system  of  pointers.  The 
array  numoff  is  kept  up  to  date  by  adding  1 to  the  row  in 
which  each  new  nonzero  term  caused  by  the  decomposition 
of  that  node  appears.  Also,  as  in  all  prior  renumbering  in 
Part  I and  Part  II.  the  array  numoff  is  altered  by  subtracting 
1 from  the  appropriate  rows  containing  the  nonzero  off- 
diagonal  terms  of  the  node  just  renumbered.  If.  by  this  sub- 
traction. an  effective  number  of  1 off-diagonal  terms  appears 
in  any  of  the  nonrenumbered  rows,  that  row  is  immediately 
renumbered  next.  After  the  bookkeeping  operations  have 
been  completed  for  renumbering  of  a row  from  Part  III, 
Part  11  is  entered  at  the  beginning.  The  search  proceeds  from 
this  point  as  if  it  were  the  first  entry  into  (I.  Fig.  1 illustrates 
the  matrix  Y at  an  intermediate  Ytage  of  renumbering.  The 
X represent  the  locations  of  all  nonzero  terms  of  an  ad- 
mittance matrix  for  a 23-node  circuit.  The  numbers  around 
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Fig.  2.  Put  I of  the  node  renumbering  algorithm. 


numbers  appear  is  the  arrangement  for  decomposition  at 
an  intermediate  stage  of  the  renumbering.  The  numbers  at 
the  far  right  are  the  current  numbers  in  the  array  numoff. 
Note  that  in  the  first  ten  rows  these  numbers  indicate  the 
number  of  nonzero  off-diagonal  terms  in  the  upper  tri- 
angular portion  of  the  array.  The  renumbering  is  complete 
to  this  tenth  row  which  lies  just  above  the  dashed  horizontal 
line.  The  two  dashed  lines  outline  a square  lower  right  comer 
array  which  is  referred  to  as  the  nonrenumbered  portion  of 
the  matrix.  The  extreme  right-hand  numbers  opposite  the 
rows  in  this  smaller  array  indicate  the  number  of  nonzero 
off-diagonal  terms  that  lie  within  this  lower  right  comer 
array.  If  any  nonrenumbered  node  was  renumbered  11,  it 
would  have  this  number  of  terms  in  the  upper  triangular 
matrix  II. 

For  the  tenth  row,  node  3 was  selected  instead  of  node  10 
because  either  one  would  cause  two  new  nonzero  terms  dur- 
ing decomposition,  but  node  3 had  more  off-diagonal  terms 
in  the  nonrenumbered  portion.  The  two  zeros  at  the  inter- 
sections of  nodes  8 and  3 indicate  the  two  locations  for 
these  new  nonzero  terms. 

After  all  of  the  nodes  are  renumbered  into  the  order  in 
which  the  nodal  equations  will  appear  in  the  admittance 
matrix,  the  stage  one  pointers  are  reorganized.  In  the  reor- 
ganization all  of  the  pointers  are  changed  to  correspond  to 
the  new  system  of  equation  numbers  and  include  all  nonzero 
terms  that  will  ultimately  be  found  in  the  upper  triangular 
matrix  V.  Terms  of  the  lower  triangular  matrix  L are  located 
by  symmetry. 

IV.  Sparse  Matrix  Solution 

The  node  voltage  vector  can  now  be  obtained  by  direct 
solution  using  Gaussian  elimination  with  the  sparse  matrix 
technique.  Because  c f the  renumbering,  the  nonzero  terms 
of  the  admittance  matrix  are  located  in  the  columnar  array 
in  the  order  in  which  decomposition  will  proceed,  by  rows 
in  the  upper  triangular  portion  and  by  columns  in  the  lower. 
This  storage  arrangement  was  selected  to  allow  for  very 
efficient  decomposition  processing  and  completely  elim- 
inates the  need  for  any  further  search  for  nonzero  terms. 


Once  all  the  values  of  L and  U are  determined  by  the 
decomposition,  the  process  of  forward  and  back  substitu- 
tion is  employed  to  obtain  the  solution  vector  x.  The  steps 
involved  are  the  following  : 

Yx  - b 

lu  - r 

Lz  = b,  i z,:  = L~'b 

Ux  » z.  i.e_  x =*  U-'L-'b. 

The  IT 1 and  U~ 1 matrices  are  never  really  determined  in 
solving  for  the  t or  t vectors  just  mentioned.  Forward 
substitution  is  used  to  find  z and  back  substitution  is  used 
to  determine  x.  A slight  modification  to  the  standard 
forward  substitution  technique  is  convenient  because  the 
nonzero  terms  of  the  L matrix  are  stored  by  column.  The 
modified  forward  substitution  process  proceeds  by  columns 
instead  of  by  rows  and  results  in  the  same  number  of  opera- 
tions as  the  row  approach.  Elements  of  the  z vector  replace 
elements  of  the  b vector  as  they  are  determined.  The  follow- 
ing is  an  example  of  the  modified  forward  substitution  tech- 
nique for  a 4 x 4 full  matrix : 

I 0 0 

I,,  1 0 

Iji  (ji  1 

Ui  Uz  Ui 

The  first  column  step  is 

!i  “ hi 

*v*  - *a  - 

*JJ>  = ^3  “ 

*<n  _ u i . 

U|-|- 

The  second  column  step  is 

*■2  zi 

w<2)  _ -t  1 1 

“3  35  -J 

_ «in 

*4  -4 


~ tsi:z 
~ 1*2*2- 
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Fig.  3.  Part  11  of  the  node  renumbering  algorithm. 


The  third  column  step  is 

*3  - *?' 

*?'  - rf'  - !«*,- 

The  fourth  column  step  is 

. 3 .131 

■*■  * ~ * 


Just  as  in  the  sparse  matrix  decomposition,  only  those 
operations  that  cause  a finite  change  in  a stored  value  in  the 
b vector  are  carried  out.  . 

The  back  substitution  used  to  determine  the  final  solution 
vector  x is  a standard  technique  proceeding  from  the  final 
nth  row  equation  backward  to  the  first  row.  The  technique 
for  a 4 x 4 full  matrix  example  is  illustrated  as : 


“tl 

“13 

“13 

“t» 

|*i' 

M 

0 

“’.2 

“S3 

“l« 

X; 
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-4 

In  general 
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Fig.  4.  Part  III  of  the  node  renumbering  algorithm. 
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Fig.  5.  Subroutine  renumber  called  from  Part  I and  Part  II  of 
the  node  renumbering  algorithm. 


Proceeding  backward  from  the  fourth  equation  we  get  the 
following  answers : 


As  in  all  previous  sparse  matrix  modifications  to  the 
standard  techniques,  only  those  multiplication  and  iddi- 
tion  operations  involving  nonzero  quantities  are  carried 
out.  The  elements  of  x replace  (he  elements  of  ; in  storage 
as  they  are  determined. 


V Summary 

The  time  savings  that  have  resulted  from  using  the  de- 
scribed sparse  matrix  techniques  in  lieu  of  a standard  lull 
matrix  algorithm  is  very  significant  for  obtaining  a fre- 
quency response  to  even  relatively  small  circuits  in  the  15- 
to  .lO-nodc  category  The  structure  of  nonzero  terms  and 
the  node  renumbering  is  determined  just  once  for  any  given 
circuit.  The  most  significant  time  savings  result  when  many 
solutions  for  cither  different  circuit  parameter  values  or 
dilferent  excitation-  frequencies  are  required.  An  ideal 
application  for  the  sparse  matrix  and  renumbering  algo- 
rithm arises  in  automated  network  design  in  which  a search 
for  optimum  circuit  parameter  values  is  carried  out  in  the 
frequency  domain  [5 1. 

The  standard  triangular  decomposition  by  Gaussian 
elimination  reauires  on  the  order  of  rr'  3 multiplication- 
addition  operations.  The  sparse  matrix  decomposition  re- 
quires some  number  proportional  to  n operations  hut  is 
very  dependent  upon  the  nonzero  structure  of  a given  cir- 
cuit. Empirical  results  based  upon  the  analysis  of  typical  20- 
to  ''O-node  integrated  circuit  amplifiers  has  indicated  this 
number  varies  from  4«  to  I6n.  but  other  circuit  examples 
that  would  fall  both  above  and  below  this  range  could  be 
found. 
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W«  TURN 

Fig.  6.  Subroutine  insert  called  from  Pan  1 1 of  the  node 
renumbering  algorithm. 

The  forward-substitution  and  back-subsutution  opera- 
tions are  also  proportional  to  n and  dependent  upon  non- 
zero structure.  Nearly  one  half  of  the  total  solution  effort 
is  expended  in  these  last  two  algorithms.  Additional  time 
savings  of  up  10  a factor  of  2 could  be  achieved  under  some 
circumstances,  if  the  node  renumbering  also  takes  into  ac- 
count the  sparse  solution  efficiency  of  the  forward  substi- 
tution and  back  substitution.  Many  times  in  a frequency 
response  only  one  or  two  node  voltages  are  desired.  If  these 
nodes  are  numbered  last,  the  back  substitution  could  be 
stopped  in  its  first  stages  after  these  desired  voltages  are 
determined.  Normally  only  one  excitation  source  exists  in 
the  circuit.  The  b vector  (current  vector)  consists  of  nearly 
ail  zero  values.  If  the  nodes  associated  with  the  nonzero  b 
values  were  to  be  numbered  next  to  last,  the  forward  sub- 
stitution could  begin  near  the  right-hand  side  of  the  L 
matrix  and  need  only  proceed  through  the  final  few  nodes. 
Unfortunately,  by  arbitrarily  numbering  certain  nodes 
last,  the  decomposition  efficiency  will  most  likely  suffer. 
However,  it  has  been  shown  experimentally  that  the  total 
solution  effort  and  the  computer  time  will  be  reduced  by 
one  third  to  one  half  in  most  cases  when  this  scheme  is  em- 
ployed. 


Fig.  7.  Example  I.  (a)  The  circuit,  tbl  The  admittance  matrix  using  the 
unctrded  node  numben.  icl  The  admittance  matrix  using  the  circled 
node  numben  determined  by  the  renumbering  algorithm. 


Appendix  I 

The  node  renumbering  algorithm  is  broken  down  into 
the  three  parts  described  in  Section  III.  Row  charts  for 
these  three  parts  are  given  in  Figs.  2— J.  The  two  subroutines 
called  from  these  How  charts  are  shown  in  Figs.  5 and  6. 
The  computer  variables  are  detined  as  follows: 


LOAD 

,V 

[ORDER  (j) 

node)/) 

ILK!  /) 

IL'Cl/) 

NUMOFFl j) 


equals  the  next  node  number  to  be  assigned 
during  renumbering;  starts  at  1 and  pro- 
ceeds sequentially  to  V. 
number  of  nodes  in  the  circuit  excluding  the 
reference  node. 

array  containing  the  original  node  numbers 
in  the  order  in  which  decomposition  will 
proceed. 

array  complementary  to  iorderI/),  i.e..  if 
IORDER(y)  = lt,  NODE! It) “ J. 
array  indicating  the  starting  location  of 
terms  in  HJC(k)  associated  with  row  j of  the 
matrix. 

array  indicating  the  column  location  of  the 
term  stored  in  /tflOO-cy)  or  the  row  location 
of  the  term  stored  in  .4(500  +j). 
array  indicating  the  effective  number  of 
nonzero  off-diagonal  terms  left  in  row  j at 
any  simulated  stage  of  the  decomposition; 
the  j is  the  ortginal  node  number. 
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Fig.  S Example  2.  lai  The  circuit,  fh'  The  ailnvHano*  rnninx  u-one  the  uncircled  nrnte  numbers  ici  The  jdmmance 
matrix  using  the  circled  node  numbers  determined  bv  the  renumbering  algorithm. 


lFILL(y)  array  indicating  the  number  ol’  new  nonzero 
terms  that  would  be  created  in  ihe  L or  I ' 
matrix  if  currently  numhered  row  > was  to  be 
renumbered  as  row  number  i.oai>. 
inserts  flag  indicating  whether  or  not  insertion  into 
the  pointer  arrays  of  a potential  new  non- 
zero term  caused  by  decomposition  should 
be  earned  out:  if  it  equals  I then  the  inser- 
tion is  done,  otherwise  it  is  not. 
all  others  dummy  integer  vanables  and  integer  arravs 

Appendix  II 

Figs.  7 and  8 are  two  circuit  examples  that  compare  the 
computational  effort  required  for  triangular  decomposition 
by  Gaussian  elimination  using  a full  matrix  algorithm.  Ihe 


sparse  mainx  algorithm,  and  the  sparse  matnx  algorithm 
with  the  node  renumbering  aigonthm.  In  noth  circuits  the 
reference  node  is  numhered  zero  and  all  terms  in  the  ad- 
mittance matrix  associated  with  the  zero  node  are  dropped. 
The  resulting  nonsingular  matrix  is  pamtioned  as 


where 

Y„  square  array  of  nodal  admittance  terms  cor- 
responding to  nodes  of  unknown  voltage 
K,  square  array  of  nodal  admittance  terms  cor- 
responding to  nodes  of  grounded  indepen- 
dent voltage  sources . 
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Y„  and  Ym  intersections  of  the  above  two  sets; 
u„  and  i.  unknown  voltage  and  independent  current 
source  vector  associated  with  the  nodes; 
u,  and  i,  voltage  source  and  current  source  vector 
associated  with  the  Y„  nodes ; 
i,  unknown  current  vector  of  currents  through 
the  grounded  independent  voltage  source. 

The  node  voltage  vector  v,  is  to  be  determined  by  Gaus- 
sian elimination  on  the  set  of  equations 

O.  * - Ym°r 

It  is  the  T_  matrix  that  is  shown  in  the  two  examples. 
Figs.  7 and  8.  The  X indicate  the  nonzero  terms  of 
before  decomposition.  The  zeros  indicate  terms  which  are 
zero  in  Y„  but  change  to  finite  values  during  the  decomposi- 
tion into  LU  form.  In  example  1,  Fig.  7,  the  standard  full 
matrix  decomposition  requires  55  divisions,  385  multipli- 
cations. and  385  additions.  By  using  the  sparse  matrix 
technique  without  renumbering,  the  decomposition  requires 
19  divisions,  41  multiplications,  and  41  additions  With 


renumbering  before  sparse  decomposition,  these  operations 
are  reduced  to  15  divisions,  25  multiplications,  and  25 
additions. 

In  example  2.  Fig.  8.  the  full  matrix  decomposition  re- 
quires 378  divisions,  6930  multiplications,  and  6930  addi- 
tions. This  compares  with  134  divisions.  746  multiplica- 
tions, and  746  additions  for  sparse  decomposition  without 
renumbering  and  63  divisions.  169  multiplications,  and  169 
additions  with  renumbering  before  sparse  decomposition. 
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USFPS  GUIDE  for  SPICE  l 

SPICF  IS  A GENERAL  PURPOSE  CIRCUIT  SIMULATION  PROGRAM  FOR  NONLINEAR  DC, 
NCNLINFAR  TRANSIENT.  AND  LINFAR  AC  ANALYSIS.  CIRCUITS  MAY  CONTAIN  RESISTORS, 
CAPACITORS,  INDUCTORS,  INDEPENDENT  VOLTAGE  ANO  CURRENT  SOURCES,  VOLTAGE 
CEPPNOFNT  CURRENT  SOURCES.  ANf  THE  FOUR  MOST  COMMCN  SEMICONDUCTOR  DEVICES: 
BJTS,  DIODES,  JFFTS  AND  MCSFSTS. 

SPICF  HAS  BUILT-IN  MODELS  FCR  Th£  SEMICONDUCTOR  DEVICES,  AND  THE  USER 
SPECIFIES  ONLY  THE  PFPTINENT  MOOEl  PARAMETER  VALUES.  TWO  MODELS  ARE  AVAILABLE 
FOR  THF  BJT.  THE  SIMPLER  MODEL  IS  BASED  ON  THF  E8ER  S-vOLL  MODEL  ANT  INCLUDES 
CHARGE  STORAGF  effects,  OHMIC  RESISTANCES,  ANO  A CURRENT  OFPFNrFN  T CUTPUT 
CONCLCTANCF . a MODEL  BASED  CN  THE  INTEGRAL  CHARGE  MODEL  OF  GUMMfl  AND  PCON  IS 
ALSO  AVAILABLE  FOR  PROBLEMS  WHICH  REQUIRE  A MORE  SOPHISTICATED  PJT  MCDFL . THF 

Otaoe  MCPFL  CAM  RE  used  FOP  either  junction  oiooes  or  shottky  barrier  oiodes. 

THE  jF^r  AND  MOSFET  MODELS  AR F BOTH  9ASED  ON  ThE  FET  MODEL  OF  SHICHMAN  AND 
HODGES. 

PROGRAM  LIMITATIONS 


2CC  NOCFS.  INCLUDING  INTF°NAL  OEVICF  NODES.  EACH  NONZERO  CHMtC  RESISTANCE 
IN  A DEVICE  WILL  GENERATE  AN  INTERNAL  NODF.  FOR  EXAMPLE,  A CIRCUIT 
WITH  35  USER  SPECIFIED  NODES  AND  10  RJTS  WITH  NCNZFRC  BASF  AND 
COLLECTOR  RESISTANCES  WILL  CCNTAIN  55  NOCES. 

50  n<=V!C.FS  IBJTS,  DIODES,  JFETS,  AND  MUSFETSI. 

?5  I N DEP c NO,f NT  VOLTAGE  OR  CUORENT  SOURCES.  ONLY  5 INDEPENDENT  SOURCES  CAN 
R F T I Mr  OEPFNDENT  FOR  TRANSIENT  ANALYSIS. 

ZDO  TOTAL  FLEMENTS,  INCLLOING  DEVICES  ANO  I NOF  °F  NOENT  SOUPCcS. 

ID  CUTPUT  VARIABLES.  AN  OUTPUT  VARIABLE  IS  E I THF  R ^A  NODF  TO  NCOF  VOLTAGE 

OR  A CURRFMT  TH=CUGH  AN  INDEPENDENT  VGLTAGF  SOUR*CF.  CLTPUT  VARIABLES 
MAY  AC  PRINTED  IN  TABULAR  FORM,  PLOTTED  AS  LINE  PRINTER  PLOTS,  OR  BOTH. 
CNLY  5 OUTPUT  VARIABLFS  CAN  BE  USED  IN  THE  AC  SMALL  SIGNAL  ANALYSIS. 

10  Sc  T S OF  MQCEL  PARAMFTFRS  FOR  DEVICFS. 
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TYPES  OF  ANALYSIS 


CC  ANALYSIS 

THE  OC  ANALYSIS  PORTION  OF  SPICE  CETFRMINES  THF  OC  OPERATING  POINT  OF  THE 
CIRCUIT  WITH  INDUCTORS  SHORTFD  AND  CAPACITORS  OPENEO.  A OC  ANALYSIS  IS 
AUTCMAT  ICALLY  PERFORMED  PRIOR  TO  A TRANSIENT  ANALYSIS  TO  DETERMINE  THE  TRANSIENT 
INITIAL  CONDITIONS.  AND  PRIOR  TO  AN  AC  SMALL  SIGNAL  ANALYSIS  TO  DETERMINE  THE 
L I NFAR I ZED.  SMALL  SIGNAL  MOOELS  FOR  NONLINEAR  DEVICES.  IF  REQUESTFC,  THE  OC 
SMALL  SIGNAL  VALUE  OF  A TRANSFFR  FUNCTION  (RATIO  CF  OUTPUT  VARIABLE  TO  INPUT 
SCLRCF),  INPUT  RESISTANCE,  AND  CUTPUT  RESISTANCE  WILL  ALSO  3E  CCMPUTED  AS  A PART 
OF  THE  SMALL  SIGNAL  OPERATING  POINT.  THE  OC  ANALYSIS  CAN  ALSO  0F  USED  TO 
GFNFRAT6  CC  TRANSFER  CURVES.  A SPECIFIED  INOEPENOENT  VOLTAGF  OP  CURRENT  SOURCE 
IS  STFPPED  OVER  A USED  SPECIFIED  RANGE  AND  THE  OC  OUTPUT  VARIABLES  ARE  STORED 
FOR  EACH  SEOUFNT  IAL  SOURCE  VALUE.  THE  OC  ANALYSIS  OPTIONS  ARE  SPECIFIED  ON  THF 
.DC  CONTROL  CARO  (PAGF  15). 

AC  SMALL  SIGNAL  ANALYSIS 

THE  AC  SMALL  SIGNAL  PORTION  CF  SPICE  COMPUTES  THE  AC  OUTPUT  VARIABLES  AS  A 
FUNCTION  OF  FREQUENCY.  THE  PROGRAM  FIRST  COMPUTES  THE  OC  OPERATING  POINT  OF  THE 
CIRCUIT  ANC  DETERMINES  LINEARIZED,  SMALL  SIGNAL  MCOELS  FCR  ALL  CF  THE  NONLINEAR 
OEVICFS  IN  THE  CIRCUIT.  THE  RESULTANT  LINEAR  CIRCUIT  IS  THFN  ANALYZED  OVER  A 
LSER  SPECIFIED  RANGE  OF  FREQUENCIES.  THE  DESIRFD  OUTPUT  OF  AN  AC  SMALL  SIGNAL 
ANALYSIS  IS  USUALLY  A TRANSFER  FUNCTION  (VOLTAGE  GAIN,  TRANS  I MPEOANCE , ETC).  IF 
THE  CIRCUIT  HAS  ONLY  ONE  AC  INPUT,  IT  IS  CONVENIENT  TO  SET  THAT  INPUT  TO  UNITY 
AND  ZERO  PHASE.  SO  THAT  OUTPUT  VARIABLES  H A V F THE  SAME  VALUE  AS  THE  TRANSFER 
FUNCTION  OF  THF  OUTPUT  VARIABLE  WITH  RFSPECT  TO  THE  INPUT. 

THE  GENERATION  OF  WHITF  NOISE  BY  RESISTORS  AND  SEMICONDUCTOR  DEVICES  CAN 
ALSO  BE  SIMULATED  WITH  THE  AC  SMALL  SIGNAL  PORTION  OF  SPICE.  EQUIVALENT  NOISF 
SOURCE  VALUES  ARE  OETERMINFO  AUTCMAT ICALLY  FROM  THE  SMALL  SIGNAL  OPERATING  POINT 
OF  THF  CIRCUIT,  ANO  THE  CONTRIBUTION  OF  EACH  NOISE  SOURCE  IS  ACCEC  AT  A GIVEN 
SUMMING  POINT.  THE  TOTAL  CUTPUT  NOISE  LEVEL  AND  THE  EQUIVALENT  INPUT  NOISE 
LEVEL  ARE  DETERMINED  AT  EACH  FREQUENCY  POINT.  THE  OUTPUT  ANO  INPUT  NOISE  LEVELS 
ARE  NORMALIZED  WITH  RESPECT  TO  THE  SQUARE  ROOT  OF  THE  NOISE  E ANC  W 1 0 TH  AND  HAVE 
THF  UMTS  VOLTS/RT  HZ  OR  AMPS/RT  HZ.  THE  OUTPUT  NOISE  ANO  EQUIVALENT  INPUT 
NOISE  CAN  BE  PRINTED  OR  PLOTTED  IN  THE  SAME  FASHION  AS  OTHER  OUTPUT  VAP  TABLES. 

THE  FPCQUFNCY  R.ANGF  ANC  THE  NOISE  ANALYSIS  OPTIONS  ARE  SPECIFIED  ON  THE 
.AC  CONTROL  CARO  (PAGE  15). 

TRANSIENT  ANALYSIS 

THE  TRANSIENT  ANALYSIS  PORTION  OF  SPICE  COMPUTES  THE  TRANSIENT  OUTPUT 
VARIABLES  AS  A FUNCTION  OF  TIMF  OVER  A USER  SPECIFIFO  TIME  INTERVAL.  THE 
INITIAL  CCNO IT  IONS  ARF  AUTOMATICALLY  DETERMINED  BY  A OC  ANALYSIS.  ALL  SOURCES 
WHICH  ARE  NOT  TIME  DEPENDENT  (FOR  EXAMPLE,  POWER  SUPPLIES)  ARE  SET  TO  THF IR  DC 
VALUE.  FOR  LARGE  SIGNAL  SINUSOIDAL  SIMULATIONS,  A FCURIER  ANALYSIS  OF  THE 
OUTPUT  WAVEFORM  CAN  BE  SPECIFIEO  TO  OBTAIN  THE  FREQUENCY  OOMAIN  FOURIER 
COEFFICIENTS.  THE  TRANSIENT  TT*E  INTERVAL  ANO  THE  FOURIER  ANALYSIS  OPTIONS 
ARE  SPECIFIED  UN  THE  .TRAN  CONTROL  CARO  (PAGE  16). 
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ANALYSIS  AT  DIFFERENT  TEMPERATURES 

ALL  INPUT  DATA  FOR  SPICE  IS  ASSUMEn  TO  HAVE  BEEN  MEASURED  AT  27  CEG  C 
(300  CFG  K).  the  SIMULATION  ALSO  ASSUMES  A NOMINAL  TEMPERATURE  OF  27  DEG  C. 

THE  CIRCUIT  CAN  dE  SIMULATED  AT  UP  TO  5 DIFFERENT  TEMPERATURES  PY  LSING  A .TE“P 
CCMRCL  CARD  I PAGE  151. 

' 

TEMPERATURE  appears  fxplicitly  IN  ThF  EXPONENTIAL  TERMS  of  the  8 JT  ano 
OIODF  MODEL  EQUATIONS.  IN  ADDITION,  SATURATION  CURRENTS  HAVE  A BUILT-IN 
TEMPERATURE  DFPFNOFNCE.  THE  TEMPERATURE  OEPENOENCE  CF  THE  SATURATION  CURRENT 
IN  ThF  RJT  MODE l S IS  DETERMINED  By: 

IS  (TFMPI  * 10  * <TEMP**3)  * EXP  (-Q  * EG  / (K  • TEMPIJ, 

WHERE  K IS  BOLTZMANS  CONSTANT,  Q is  the  ELECTRONIC  CHARGE.  10  IS  A CONSTANT,  ANO 
FG  IS  THE  ENERGY  GAP  WHICH  IS  A MODEL  PARAMETER.  THE  TEMPERATURE  OEPFNCENCF  OF 
THf  SATURATION  CURRENT  IN  THE  JUNCTION  CIOOE  MOOEL  IS  OETERMI NFC  BY: 

IS  (TEMP)  » 10  * ( TEMP** ( 3 /N  II  * EXP  (-Q  * EG  / (K  * TFMF)), 

• 

WFFRE  N IS  the  EMISSION  COFFFICIENT,  WHICH  IS  A MCDEL  PARAMETER,  AN'C  THE  OTHER 
SYMBCLS  HAVE  THE  SAME  MEANING  AS  ABOVE.  FOR  SHOTTKY  BARR  I FR  CICTFS,  THE 
TEMPERATURf  OEPFNCENCF  OF  THE  SATURATION  CURRENT  IS  DETERMINED  BY: 

IS  (TEMP)  a io  * ( TFMP** ( 2/N ) ) * EXP  (-Q  * EG  / (K  * TEMP)). 


CONVERGENCE 


BOTH  DC  ANO  TRANS  I FNT  SOLUTIONS  ARE  OBTAINED  BY  AN  ITERATIVE  PROCESS  WHICH 
IS  TERMINATED  WHEN  THE  NOOE  VOLTAGES  CONVERGE  TO  WITHIN  A TOLFRANCF  OF  0.1 
PERCENT  OP  50  MICROVOLTS.  WHICHEVER  IS  LARGER.  ALTHOUGH  THE  PARTICULAR 
ALGORITHM  USED  IS  SPICE  HAS  B6FN  FOUND  TO  BE  VERY  RELIABLE.  IN  SCMF  CASES  IT 
WILL  FAIL  TO  CONVERGE  TO  A SOLUTION.  WHEN  THIS  HAPPENS,  THE  PROGRAM  WILL  PRINT 
OUT  THF  LAST  NOOE  VOLTAGES  ANC  TERMINATE  THE  JOB.  THE  NOOE  VOLTAGES  THAT  ARE 
PRINTED  APE  NOT  NECESSARILY  CORRECT  OR  EVFN  CLOSE  TO  THE  CORRECT  SOLUTION. 

FAILURE  TO  CONVERGF  IN  THE  OC  ANALYSIS  IS  USUALLY  DUE  TC  AN  FRROR  IN 
SPECIFYING  CIRCUIT  CONNECTIONS,  ELEMENT  VALUES,  OR  MODEL  PARAMETER  VALUES. 
REGENERATIVE  SwITCHING  CIRCUITS  OR  CIRCUITS  WITH  POSITIVE  FEECPACK  PROBABLY  WILL 
NCT  CONVERGE  IN  THE  OC  ANALYSIS.  FAILURE  TO  CONVERGE  IN  THE  TRANSIENT  ANALYSIS 
CAN  ALSO  BE  OUS  TO  A TIME  STEP  WHICH  IS  TOO  LARGE.  SPICE  PRESENTLY  DOES  NOT 
have  AN  AUTOMATIC  time  STEP  CONTROL,  ANO  SIGNIFICANT  ERROR  ANO /CP  NCNCONVERGeNCE 
CAN  RFSULT  IE  THE  TIME  STEP  IS  LARGE  COMPARED  TO  THE  CIRCUIT  TIME  CONSTANTS. 
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INFLT  FORMA  T 


THE  INPUT  FORMAT  FOR  SPICE  IS  OF  The  FREE  FORMAT  TYPE.  FIFCCS  ON  A CArC 
ARE  SEPARATED  -i y QNF  OR  MORE  BLANKS,  A COMMA,  OR  AN  EQUAL  <«)  SIGN.  SPACES 
PRECEOING  OR  FOLLOWING  A COMMA  OR  EQUAL  SIGN  ARE  IGNCREO.  A CAPO  may  BE 

contjnlfo  onto  the  following  caho  by  punching  a ♦ before  the  first  field  on  thf 

CONTINUATION  CARO. 

A NAME  F I c LO  ML S T BEGIN  WITH  A LETTER  (A  THRU  Z)  ANC  CANNOT  CONTAIN 

commas  or  blanks,  only  the  first  seven  characters  cf  the  name  are  lsfd. 

A NUMRFR  FIELO  Mar  BE  AN  INTEGER  FIFLO  (l 2.-44),  A FLOATINC  POINT  FIELD 
(3.LBL59),  CITHER  AN  INTEGER  OR  A FLOATING  POINT  NUMBER  FOLLOWED  BY  AN  INTEGER 
FXFCNFNT  (1F-14,  2.65E3I,  OR  EITHER  AN  INTEGER  OR  A FLOATING  POINT  NUMBER 
FCLLCWFO  BY  ONE  OF  THE  FOLLOWING  SCALE  FACTORS: 

G I.CE9 

MEG  1.0F6 

K 1 . OF  3 

M l.OE-3 

U l.OE-6 

N l.OE-9 

P 1.0E-12 

LETTERS  IMMEDIATELY  FOLLOWING  A NUMBER  THAT  ARE  NOT  SCALE  FACTORS  ARE  IGNORED, 
ANO  LETTERS  [MMECIA T ELY  FOLLOWING  A SCALE  FACTOR  ARE  IGNORED.  HENCE , 10,  tOV, 
10VCLTS,  ANO  10HZ  ALL  REPRESENT  THE  SAME  NUMBER,  ANO  M,  MA,  MSEC,  ANO  MMHOS  ALL 
REPPESENT  THE  SAME  SCALE  FACTOR.  NOTE  THAT  1000,  10C0.0,  1000HZ,  IF3,  1.0E3, 
IKHZ,  ANO  IK  ALL  REPRESENT  THE  SAME  NUMBER. 


CIRCLIT  DESCRIPTION 


THE  CIRCUIT  TO  BE  ANALYZED  IS  OESCRIBEC  TO  SPICE  BV  A SET  CF  ELEMENT 
CARCS.  WHICH  DEFINE  THE  CIRCUIT  TOPOLOGY  ANO  ELEMENT  VALUES.  ANC  A SET  OF 
CONTROL  CAROS,  WHICH  DEFINE  THE  MCOEL  PARAMETERS  ANO  THE  RUN  CONTROLS.  THE 
FIRST  CARO  IN  THE  INPUT  DECK  MUST  BE  A TITLE  CARD,  ANC  THE  LAST  CARC  MUST  BE  A 
•ENC  CARO.  THE  OROER  OF  THE  REMAINING  CAROS  IS  ARBITRARY. 

NOOE  NUMBERS  MUST  BE  INTEGERS.  THE  DATUM  NOOE  HUST  BE  NUMBERED  0 I ZERO). 
NODES  NEEO  NOT  BE  NUMBEREO  SEQUENTIALLY.  THE  CIRCUIT  CANNOT  CONTAIN  A LOOP  OF 
VOLTAGE  SOURCES  ANO/OR  INDUCTORS  AND  CANNOT  CONTAIN  A CUTSET  OF  CURRENT  SOURCES 
ANC/OR  CAPACITORS.  EACH  NOOE  IN  THF  CIRCUIT,  INCLUDING  THE  DATUM  NOOE,  MUST 
HA VF  AT  LEAST  TWO  CONNECTIONS. 


■ A 


*****  RESISTORS,  capacitors,  inductors 

general  FDPM  RXXXXXX  NL  N2  value 

CXXXXXX  Nl  N2  value 
LXXXXXX  Nl  N 2 VALUE 

EXAMPLES  R 13  12  17  IK 

CGOOO  13  0 10P 
LL INKS  42  69  IU 

Nl  ANO  N2  ARE  THE  TWO  ELEMENT  NODES.  THE  ORCER  OF  THE  NODES  FOR  THESE 
ELEMENTS  IS  UNIMPORTANT.  VALUE  IS  THE  RESISTANCE  IOHMSI,  THE  CAPACITANCE 
(FARACSI,  ANO  THE  IN'OUCTANCE  (HENRIES).  RESPECT  I VEL  Y . THIS  VALUE  CANNOT  BE 
NEGATIVE  OP  ZERO. 


*****  VOLTAGE  CONTROLLED  CURRFNT  SOURCES 

GENEFAL  FCRM  I XXX XXX  V N+  N-  NC*  NC-  VALUE  DELAY 

EXAMPLES  ISORS  V 13  12  14  12  l.OM 

IGM  V 1 20  4 20  -2.0M  3.0NS 

THE  LETTER  V MLST  BE  IN  THE  FIELD  FOLLOWING  THE  ELEMENT  NAME.  N*  AND  N- 
ARE  THE  POSITIVE  ANO  NEGATIVE  NODES,  RESPECT IVELY . CURRFNT  FLOWS  FROM  THE 
POSITIVE  NODE.  THRU  THE  SOURCE.  TC  THE  NEGATIVE  NOOE.  NC>  ANO  NC-  ARE  THE 
POSITIVE  ANO  NEGATIVE  CONTROLLING  NOCES,  RESPECT  I VEL Y.  VALUE  IS  THE 
TR ANSCONOUCTANCE  (MHOS). 


IN  THE  AC  ANALYSIS  THE  T PANS CONDUCTANCE  CAN  BE  MOOIFIEC  BY  AN  OPTIONAL 
OELAY  (LINEAR  PHASE)  OPERATOR.  ThE  Oei‘AY  (SECONOS)  IS  APPENDED  AFTER  THE  VALUE. 
IF  A DELAY,  TO,  IS  INCLUDED,  THE  COMPLEX,  FREQUENCY  CEPENOENT  VALUE  OF 
TRANSCCNOUCTANCF  IS  OETERMINEC  BY: 


GM  * VALUE  * EXP  (-J  * 6.28318  * FREQ  * TO) 


THE  OELAY  IS  IGNORED  IN  THF  DC  ANO  TRANSIENT  ANALYSES. 
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**•*•  INOFPENOFNT  SOURCES 

CFNERAL  ETRM  VXXXXXX  N + N-  CC  DCVAl  AC  ACVAL  PHASE 

IXXXXXX  N*  N-  CC  OCVAC  AC  ACVAL  PHASE 

FXAMPLES  VCC  10  0 DC  6 

IZENER  13  15  DC  600MA 
VIN  13  2 DC  0.001  AC  l 
1 1 N 21  23  AC  0.333  A5.0 
VMEAS  12  9 

N>  IS  THF  POSITIVE  NOPE  ANO  N—  IS  THE  NEGATIVE  NODE.  NOTE  THAT  VOLTAGF 
SOURCES  NFFC  NOT  RE  GROUNDED.  CURRENT  FLOWS  FROM  THE  POSITIVE  NODE,  THRU  THE 
SOURCE,  TO  THE  NEGATIVE  NOOE. 

DCVAL  IS  THF  DC  VALUE  OF  THE  SOURCE.  THE  SOURCE  IS  SET  TO  THIS  VALUE  FOR 
DC  ANALYSIS  ANO.  IF  NO  TIME  DEPENDENCE  IS  ATTACHED.  IN  THE  TRANSIENT  ANALYSIS. 

IF  THF  DC  SOURCE  VALUE  IS  ZERO,  THE  LETTERS  DC  ANO  THE  OC  VALUE  CAN  BE  OMITTED. 

ACVAL  IS  THE  AC  VALUE  ANC  PHASE  IS  THE  AC  PHASE.  THE  SOURCE  IS  SET  TO  THIS 
VALUE  IS  THE  AC  ANALYSIS.  THE  ARBITRARY  PHASE  FACTOR  CAN  BE  OMITTFO.  IF  THE 
SOURCE  IS  NOT  AN  AC  SMALL  SIGNAL  INPUT,  THE  LETTERS  AC  ANO  THE  AC  VALUES  ARF 
CMITTEC. 


A SOURCE  MAY  RE  GIVEN  A TIME  DEPENDENCE  FOR  THE  TRANSIENT  ANALYSIS  BY 
APPFNDING  CNF  OF  THE  THRFE  PREDEFINED  FUNCTIONS:  PULSE,  EXPONENTIAL,  ANO 
SINLSOICAL.  IF  PARAMETERS  CTFER  THAN  SOURCE  VALUES  ARE  OMITTED  OR  SET  TC  ZERO, 
THF  DEFAULT  VALUES  SHOWN  WILL  BE  ASSUMED.  T STEP  IS  THE  PRINTING  INCREMENT  ITIME 
STEPI,  ANO  T STOP  IS  THE  FINAL  TIME  (PAGE  15). 


1.  PLLSE  PULSE  VI 


EXAMPLE 

VIN 

3 0 PULSE  -l 

PARAMETERS  ANO  OEFAULT 

VALUES 

VI 

INITIAL  VALUE 

— 

V2 

PULSED  VALUE 

— 

TO 

CELAY  time 

T STEP 

TR 

RISE  TIME 

TSTEP 

TF 

FALL  TIME 

TSTEP 

PW 

WIDTH 

T STOP 

PER 

PERIOD 

TSTOP 

A SINGLE  PULSE  IS 

DESCR I BEO  BY 

TIME 

VALU 

E 

0 

VI 

TO 

VI 

TDtTR 

V2 

TOMR^PW 

V2 

to*tr*pw*tf  VI 

TS  TOP 

VI 

V2  TO  TR  TF  PW  PER 
1 2NS  2NS  2NS  50NS  100NS 


THE  FOLLOWING  PIECEWISE  LINEAR  TABLE. 
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2.  EXPONENTIAL  FXP  vt  V2  TDl  TAUl  TD2  TAU2 

FXANPL6  VIN  1 0 EXP  -A  -1  2NS  30NS  60N$  AONS 


PARAMETERS  ANO  DEFAULT  VALUES 

VI  INITIAL  VALUE  

V2  PULSED  VALUE  

TO l RISE  DELAY  TIME  TSTEP 

TAUl  RISE  TINE  CONSTANT  TSTEP 

T02  FALL  DELAY  TIME  TSTEP 

TAU2  FALL  TINE  CONSTANT  TSTEP 

TINE  VALUE 


0 TO  TO  1 VI 

101  TO  TD2  Vl*!V2-Vl M1-EXP(-(T-T01I/TAUI)  I 

T02  TO  TSTOP  V l + I V2-V U I l-FXPI - ( T-TO l I /TAU L I !♦< V l-V 2 )( l-EXPI -( T- T02 > /TAU2 I I 

3.  SINLSOIOAL  SIN  VO  VA  FREQ  TC  THETA 

EXAMPLE  VIN  3 0 SIN  0 l 100MEG  INS  1E10 

PARAMETERS  AN C OE FAULT  VALUES 

VO  OFFSET 

VA  AMPLITUOE 

FREQ  FREQUENCY  (IN  HZ  I 

TO  DELAY 

THETA  OANPING  FACTOR 

TINE  VALUE 

0 TO  TO  VO 

TD  TO  TSTOP  VO  *■  VA  * EXPl-(T-TO)  * THETA)  * SINE  (6.28318  * FREQ  * T) 


SOURCES  MAY  BE  GIVEN  ANY  COMBINATION  OF  VALUES  (OC,  AC,  CP  TRANSIFNT),  ANO 
THESE  VALUFS  MAY  BE  SPECIFIED  IN  ANY  ORDER  AS  LONG  AS  THEY  FOLLOW  THE  PROPER 
KEY  WGRO. 

EXAMPLES  VIN  13  12  SIN  0 1 10MEG  OC  0.1  AC  l *5 

IZ  ID  0 OC  0 PULSE  0 1 AC  0.5 

VEQ  12  0 OC  0.5  EXP  0.5  0.9  IONS  AONS  70NS  AONS  AC  I 


l/TSTOP 

TSTEP 

0 
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*****  BIPOLAR  JUNCTION  TRANSISTORS 


GENERAL  FORM 


0 XXX XXX  NC  N8  NE  MNAME  AREA 


EXAMPLE 


CAMP 3 3 7 9 l MODI  2.0 


NC  IS  the  collector  node,  nb  IS  The  base  node,  ne  is  the  emitter,  noce, 

MNAME  IS  the  MOOEL  NAME  (PACF  9 ) ANO  ARFA  IS  THE  AREA  FACTOR.  THE  AREA  FACTOR 
IS  EQUIVALENT  TO  THE  NUMBER  OF  PARALLEL  OEVICES.  AN  AREA  FACTOR  OF  2.0 
I “PL I E S THAT  TWO  TRANSISTORS  OF  THE  SAME  MOOEL  ARE  CONNECT EC  IN  PARALLEL.  IF 
THF  AREA  is  OMITTED,  AN  AREA  FACTOR  OF  1.0  IS  ASSUMEO. 

*****  JUNCTION  OIOCES 


GFNERAL  FORM 


OXXXXXX  N*  N-  MNAME  AREA 


EXAMPLE 


CBR10CE  8 10  0I00E1 


N*  IS  THE  POSITIVE  NOOE,  N-  IS  THF  NEGATIVE  NOOE.  MNAME  IS  THE  MOOEL  NAME 
(PACE  SI.  ANO  AREA  IS  THE  AREA  FACTOR  (SEE  8JTS . ABOVE). 


*****  JUNCTION  FIELD  EFFECT  TRANSISTORS 


GENERAL  FORM 


JXXXXXX  NO  NG  NS  MNAME  AREA 


EXAMPLE 


J1  7 2 3 JMl 


NO  IS  THF  CRAIN  NOOE.  NG  IS  THE  GATE  NODE.  NS  IS  THE  SOURCE  NOCE,  MNAME  IS 
THE  MOOEL  NAMEIPAGE  9),  ANO  AREA  IS  THE  AREA  FACTOR  (SEE  8JTS,  ABOVE). 


MO  SF  FT  S 


GFNERAL  FCRM 


MXXXXXX  NO  NG  NS  N8  MNAME  AREA 


EXAMPLF 


M3 1G  2 3 4 7 MLONG 


NO  IS  THE  CRAIN  NOOE,  NG  IS  THE  GATE  NODE,  NS  IS  THE  SOURCE  NCCE,  N8  IS  THE 
8LLK  (SUBSTRATE)  NOOE,  MNAME  IS  THE  MOCEL  NAME  (PAGE  9),  ANO  AREA  IS  THE  AREA 
FACTOR  (SEE  BJTS,  ABOVE). 
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*****  .MODEL  CARO 

GENERAL  FORM  .MODEL  MNAME  TYPE  PNAM€l*PVALl  PN4ME2=PVAL2  ... 

EXAMPLE  .MODEL  MODI  NPN  BF*50  I S=  IE—  13  VA=*50 

THE  .MODEL  CARC  SPECIFIES  A SET  OF  MODEL  PARAMETERS  THAT  WILL  BE  USED  BY 
CNE  CR  MORE  DEVICES.  MNAMF  IS  THE  MOOEL  NAME.  ANO  TYPF  IS  CNE  CF  THF  FOLLOWING 
TEN  TYPES: 

NPN  NFN  EBERS-MOLL  8JT  MODEL 

PNP  PNP  EBERS-MOLL  BJT  MODEL 

NGP  NPN  GUMMEL-POON  BJT  MOOEL 

PGP  PNP  GUMMEL-PCCN  BJT  *OOEL 

D JUNCTION  C I ODE  MODEL 

SBO  SHOTTKY  BARRIER  OIOCE  MODEL 

NJF  N CHANNEL  JFET  MOCEL 

PJF  0 CHANNEL  JFET  MCOEL 

NMO  N CHANNEL  MOSF FT  MOOEL 

PMO  P CHANNEL  MOSFET  MQCEL 

PARAMETER  VALUES  ARE  DEFINED  BY  APPENDING  THE  PARAMETER  NAME,  AS  GIVEN 
BECCW  FOR  EACH  MODEL  TYPE,  FOLLOWED  BY  AN  ECUAL  SIGN  ANO  THE  PARAMETER  VALUE. 
MOOEL  PARAMETERS  THAT  ARE  NOT  GIVEN  A VALUE  ARE  ASSIGNED  THE  OEFAULT  VALUE 
GIVEN  BELOW  FOR  EACH  MOCEL  TYPE. 

MODEL  PARAMETER  VALUES  CAN  ALSO  BE  SPECIFIED  AS  A STRING  CF  NUMBERS  IN  THE 
ORCER  GIVEN  BELOW  FOR  EACH  MODEL  TYPE.  THE  FOLLOWING  MOOEL  SPECI F TCAT TON  IS 
EQUIVALENT  TO  THE  PREVIOUS  MOCEL  CARO  EXAMPLE: 

EXAMPLE  .MODEL  MODI  NPN  50. IE-1 3 ...  50 

CIOOE  MOOELS  1 80 TH  JUNCTION  ANO  SBO) 


THE  CNLY  DIFFERENCE  BETWEEN  THE  JUNCTION  DICOE  MODEL  AND  THE  SHOTTKY 
BARRIER  OIODC  MOOEL  IS  THE  TEMPERATURE  DEPENDENCE  OF  SATURATION  CURRENT  (SEE 
PAGE  3).  THE  DC  CHARACTERISTICS  OF  THF  OIOOE  ARE  DETERMINED  BY  THE  PARAMETERS 
IS  ANC  N.  AN  OHMIC  RESISTANCE,  RS,  IS  INCLUOEO.  CHARGE  STCRAGF  EFFECTS  ARE 
MOOELEO  BY  A TRANSIT  TIME,  TT , ANO  A NONLINEAR  DEPLETION  LAYER  CAPACITANCE  WHICH 
VARIES  AS  THE  -1/2  POWER  OF  JLNCTtON  VCLTAGE  ANO  IS  CEFINED  BY  THE  PAPAmETFRS 
CJC  AND  PHI.  THE  ENERGY  GAP,  FG,  AFFECTS  ONLY  THE  TEMPERATURE  CEPFNOENCE  OF  THE 

saturation  current  isee  page  3). 


NAME 

PARAMETER 

OEFAULT 

TYPICAL 

1 

RS 

OHMIC  RESISTANCE 

0 

10 

2 

TT 

TRANSIT  TIME 

0 

C.1NS 

3 

CJO 

ZERO  BIAS  JUNCTION  CAPACITANCE 

0 

2PF 

6 

IS 

SATURATION  CURRENT 

1.06-14 

l .OF- 16 

5 

N 

EMISSION  COEFFICIENT 

l 

1.0 

6 

PHI 

JUNCTION  POTENTIAL 

1 

0.6 

7 

EG 

ENERGY  GAP 

1. 11  SI 

1.11  FDR 

0.69  SBC  0.69  FDR  SBD 
0.67  FOR  GE 


■ I 
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FBERS-MOLL  BJT  MQOELS  (807H  NPN  AND  PNP) 

THE  FBERS-MCLL  BJT  MODEL  USES  THE  OC  EBERS-MGLL  MGPEL  AS  A BASIS.  THE  CC 
CHARACTERISTICS  OF  THE  DEVICE  ARE  DETERMINED  BY  THE  PARAMETFP  S B F ANO  eR,  THF 
FORWARC  ANO  RFVFRSE  CURRENT  GAINS,  V A,  WHICH  DETERMINES  THE  OUTPUT  CONDUCTANCE, 

anc  the  saturaticn  current,  rs.  three  ohmic  resistances,  rb,  pc,  and  re,  have 

BFEN  INCLUDED.  BASE  CHARGE  STORAGE  IS  MOOElED  BY  FORWARD  AND  REVERSE  TRANSIT 
TIMES,  TF  ANO  TR,  ANO  NCNLINEAR  DEPLETION  LAYER  CAPACITANCES  WHICH  VARY  AS  THE 
-1/2  POWER  OF  JUNCTICN  VOLTAGE  ANC  ARE  CEFINED  BY  TEE  PARAMETERS  CJ f , PE,  CJC, 
ANO  PC.  A CONSTANT  COLLECTOR-SUBSTRATE  CAPACITANCE,  CCS,  IS  ALSO  INCLUOFO.  THE 
ENERGY  GAP,  EG.  AFFECTS  ONLY  THE  TEMPERATURE  DEPENDENCE  OF  THE  SATURATION 


CURRENT 

I SEE  PAGE 

3). 

NAME 

PARAMETER 

DEFAULT 

TYPICAL 

l 

8F 

FORWARD  BETA 

100 

100 

2 

8R 

REVERSE  BETA 

l 

0.1 

3 

RB 

BASE  OHMIC  RESISTANCE 

0 

100 

4 

RC 

COLLECTOR  OHMIC  RESISTANCE 

0 

10 

5 

RE 

EMITTER  OHMIC  RESISTANCE 

0 

l 

6 

CCS 

COLLECTOR-SUBSTRATE  CAPACITANCE 

0 

2PF 

7 

TF 

FORWARD  TRANSIT  TIME 

0 

C.1NS 

8 

TR 

REVERSE  TRANSIT  TIME 

0 

IONS 

9 

CJF 

ZERO  BIAS  B-E  JUNCTION  CAPACITANCE 

0 

2PF 

10 

CJC 

ZERU  BIAS  B-C  JUNCTION  CAPACITANCE 

0 

1PF 

li 

IS 

SATURATION  CURRENT 

l.OE-14 

l .OE-14 

12 

PE 

B-E  JUNCTICN  POTENTIAL 

l 

0.7 

13 

PC 

B-C  JUNCTICN  POTENTIAL 

1 

0.5 

14 

VA 

EARLY  VOLTAGE 

INFINITE 

50 

15 

EG 

ENERGY  GAP 

1.11 

i.n  for 

0.67  FOR  GE 


GUMMEL-PCCN  eJT  MODELS  (BOTH  NPN  ANO  PNP I 

THE  INTEGRAL  CHARGE  MOCEL  OF  GUMMEL  ANO  POON  IS  A MORE  COMPLICATED  ANO  MORE 
COMPLETE  BJT  MOOEL  FOR  PROBLEMS  WHICH  REQUIRE  ACCURATE  BJT  MOOELS.  THE  OC  MOOFi. 
IS  OEFINEO  BY  THE  PARAMETERS  BFM , C2,  IK,  ANO  N6.  WHICH  DETERMINE  THE  FORWARC 
CURRENT  GAIN  CHARACTERISTICS,  8RM , C4,  IKR,  ANO  NC,  WHICH  DETERMINE  THE  REVERSE 
CURRFNT  GAIN  CHARACTERISTICS,  VA  AND  VB,  WHICH  DETERMINE  THE  OUTPUT  CONDUCTANCE 
FCR  FORWARC  ANO  REVERSE  REGIONS,  AND  THE  SATURATION  CURRENT,  IS.  THREE 
OHMIC  RESISTANCES,  R8,  RC,  ANC  RE,  ARE  INCLUDED.  BASE  CHARGE  STORAGE  IS  MOCEIEO 
BY  FORWARC  ANO  REVERSE  TRANSIT  TIMES,  TF  ANO  TR,  AND  NONLINEAR  CEPLETICN  LAYER 
CAPACITANCES  WHICH  ARE  CETERM  INFO  BY  CJE,  PE,  ANO  ME  FOR  THE  B-t  JUNCTION,  ANO 
CJC,  PC,  ANO  NC  FOR  THE  B-C  JUNCTION.  A CONSTANT  COL LECTOR- SUBSTRATE 
CAP  AC  IT  ANCF , CCS,  IS  ALSO  INCLUOEC.  THE  ENERGY  GAP,  EG,  IS  INCLUCEC  AS  IN  THE 
StEELER  BJT  MCOFL. 


pr 

PAGE  11 

NAME 

parameter 

DEFAULT 

TYPICAL 

i 

PFM 

SAT  CURRENT/IDEAL  ft-E  SAT  CURRENT 

100 

100 

2 

BRM 

SAT  CURRENT/IDEAL  8-C  SAT  CURRENT 

1 

0.  1 

3 

RB 

BASE  OHMIC  RESISTANCE 

0 

100 

4 

RC 

COLLECTOR  OHMIC  RESISTANCE 

0 

10 

5 

RE 

emitter  ohmic  resistance 

0 

1 

6 

CCS 

COLLECTOR— SU8S TR  ATE  CAPACITANCE 

0 

2PF 

7 

TF 

FORwAPO  TRANSIT  T IMF 

0 

C.  INS 

8 

TP 

REVERSE  TRANSIT  TIME 

0 

IONS 

8 

C J F' 

ZERO  BIAS  B-E  JUNCTION  CAPACITANCE 

0 

2PF 

10 

CJC 

ZERO  BIAS  e-C  JUNCTION  CAPACITANCE 

0 

IPF 

n 

IS 

SATURATION  CURRENT 

l.OE-14 

1 .OF—  14 

12 

VA 

FORWARO  EARLY  VOLTAGE 

INFINITE 

50 

13 

VB 

REVERSE  EARLY  VOLTAGE 

INF  INITE 

50 

14 

C2 

NONIOEAL  B-E  SAT  CURRENT/SAT  CURRENT 

0 

10. 0 

15 

IK 

FORWARO  KNFE  CURRENT 

INF INITE 

10MA 

i 

16 

NE 

e-e  emission  coefficient 

2.0 

1.5 

17 

C4 

NCNI DEAL  B-C  SAT  CURRENT/SAT  CURRENT 

0 

1.0 

18 

Ik  r 

REVERSE  KNEE  CURRENT 

INFINITE 

100MA 

19 

NC 

B-C  EMISSION  COEFFICIENT 

2.0 

1.5 

20 

P F 

B-E  JUNCTION  POTENTIAL 

1 .0 

0.7 

21 

ME 

B-E  grading  coefficient 

0.5 

0.33 

22 

PC 

B-C  JUNCTION  POTENTIAL 

1.0 

0.5 

23 

MC 

8-C  GRAOING  COEFFICIENT 

0.5 

0.33 

24 

EG 

ENERGY  GAP 

1.11 

1.11  FOR  SI 

0.67  FOR  GE 

— 

JFFT  MODELS  (BOTH  N ANO  P CHANNEL) 

U'  • THE  JF  ET  MOOEL  IS  DERIVED  FROM  THE  FET  MODEL  OF  SHICHMAN  ANO 

HCDGES.  THE 

ET  CC  CHARACTERISTICS  ARE  0EFIN60  BY  THE  PARAMETERS  VTO  ANO 

beta,  WHICH  DETERMINE 

f ' THE  VARIATION  OF 

ORAIN  CURRENT  WITH  GATE  VOLTAGE,  LAMBOA 

, WHICH  DETERMINES  THE 

output 

CONDUCTANCE,  ANO  IS,  THE  SATURATION  CURRENT  OF  THE  TWO  GAT F 

JUNCTIONS.  TWO 

OHMIC  RESISTANCES,  RO  ANO  RS,  ARE  INCLUOEO.  CHARGE  STORAGE  IS 

j 

MODELED  BY  NONLINEAR  DEPLETION  LAVER  CAPACITANCES  FOR  BOTH  GATE  JUNCTIONS  WHICH 

VARY  AS  THE  -1/2 

POWER  OF  JUNCTION  VOLTAGE  ANO  ARE  DEFINED  BY  THE 

parameters 

CGS.  CGD,  ANO  P0 

• 

NAME 

parameter 

OEFAULT 

TYPICAL 

1 

VTO 

THRESHCLO  VCLTAGE 

-2.0 

-2.0 

2 

BETA 

transconouctance  parameter 

l.OE-4 

l.OE-4 

3 

LAM BOA 

CHANNEL  LENGTH  MODULATION  PARAMETER 

0 

0.01 

4 

RO 

CRAIN  OHMIC  RESISTANCE 

0 

100 

5 

RS 

SOURCE  OHMIC  RESISTANCE 

C 

100 

6 

CGS 

ZERO  BIAS  G-S  JUNCTION  CAPACITANCE 

0 

2PF 

7 

CGD 

ZERO  BIAS  G-D  JUNCTION  CAPACITANCE 

0 

2PF 

8 

PB 

GATE  JUNCTION  POTENTIAL 

1 

0.6 

9 

IS 

GATE  JUNCTION  SATURATION  CURRENT 

l.OE-14 

l.OE-14 

I 

PAGE  12 


MCSFET  MODELS  (BOTH  N ANO  P CHANNELS! 


& 


THE  MOSFET  MODEL  IS  ALSO  CERIVED  FROM  The  FET  MCCEL  OF  SHICHMAN  ANC 
HDOGFS.  THF  nc  CHARACTERISTICS  OF  THE  MOSFET  ARE  DEFINED  BY  THE  PARAMETERS 
VTC,  BFTA,  ANO  LAM6CA,  WHICH  ARE  IDENTICAL  TO  THE  PARAMETERS  FOR  THE  JFET « PHI 
ANC  CA*MA,  WHICH  DETERMINE  THE  VARIATION  OF  THRESHOLD  VOLTAGE  WITH  SUBSTRATE 
VOLTAGE.  AND  IS,  THE  SATURATION  CURRENT  OF  THE  TWO  SUBSTRATE  JUNCTIONS.  CHARGE 
STOAACF  tS  «OOELEO  BY  THREE  CONSTANT  CAPACITORS.  CGS , CGC , ANO  CCe,  ANO 
NONLINEAR  DEPLETION  LAYER  CAPACITANCES  FOR  BOTH  SUBSTRATE  JUNCTIONS  WHICH  VARY 
AS  THE  -1/2  PGwER  OF  JUNCTION  VOLTAGE  ANO  ARE  OETERMINEO  BY  THE  PARAMETERS  C8D, 
CBS,  ANC  PB. 


NAME 

P ARAMtTFR 

OEFAULT 

TYP [CAL 

1 

VTO 

THRESHOLD  VOLTAGE 

2.0 

2.0 

2 

PHI 

SURFACE  POTENTIAL 

0.5 

0.5 

3 

BETA 

TRANSCONCUCTANCE  PARAMETER 

l.OE-4 

1 . OE— 4 

4 

GAMMA 

BULK  THRESHOLD  PARAMETER 

0 

0.5 

5 

LAM6C4 

CHANNEL  LENGTH  MODULATION  PARAMETER 

0 

0.01 

6 

RD 

DRA IN  OHMIC  RESISTANCE 

0 

ICO 

7 

RS 

SOURCE  OHMIC  RESISTANCE 

0 

100 

8 

CGS 

GATE-SOURCE  CAPACITANCE 

0 

IPF 

9 

cr.n 

GATE-ORAIN  CAPACITANCE 

0 

1PF 

10 

CG8 

CATE-BULK  CAPACITANCE 

0 

IPF 

11 

CBO 

zero  eiA$  e-c  junction  capacitance 

0 

IPF 

12 

CBS 

ZERO  BIAS  e-s  junction  capacitance 

0 

IPF 

13 

PB 

BULK  JUNCTION  POTENTIAL 

1 

0.6 

14 

IS 

BULK  JUNCTION  SATURATION  CURRENT 

1.0E-14 

1.0E-14 

PAGE  13 


CONTROL  CARDS 

*****  TITLE  CARD 

EXAMPLE  OP  AMP  CIRCUIT  JOE  J STUDENT  EECS  241 

THIS  CARO  MUST  BE  THE  FIRST  CARD  IN  THE  INPUT  CECK.  ITS  CONTENTS  ARE 
PR  I NT  FC  VERBATIM  AS  THE  HEADING  FOR  EACH  SECTION  OF  OUTPUT. 

*****  .END  CARC 
GENEFAL  FORM  .END 

THIS  CARD  MtST  ALWAYS  BE  THE  LAST  CARO  IN  THE  INPUT  OECK.  NOTE  THAT  THE 
P FR  IOO  IS  AN  INTEGRAL  PART  OF  THE  NAME.  IF  THE  .ENO  CARD  IS  OMITTFC,  THE  NEXT 
JOe  WILL  8E  READ  IN  AS  PART  OF  THE  JOB  MISSING  THE  .ENO  CARO,  AND  NEITHER  JOB 
WILL  BP  RUN  SUCCESSFULLY. 

*****  COMM  PM T CARD 

GENERAL  FORM  * ANY  COMMENTS 

EXAMPLE  * RF  3 IK  GAIN  SHOULO  BE  100 

THIS  CARD  IS  PRINTED  OUT  IN  THE  INPUT  LISTING  BUT  IS  OTHERWISE  IGNCRFO . 

*****  NO  PRINT  CARO 
GENERAL  FORM  .NP 

THIS  CARO  SUPPRESSES  THF  SUMMARY  OF  INPUT  DATA  THAT  IS  NORMALLY  PRINTED 
AFTER  READING  THE  INPUT  DECK.  IT  OOES  NOT  SUPPRESS  THE  LISTING  OF  THE  INPUT 
OECK  OR  ANY  ERROR  MESSAGES  THAT  MAY  OCCUR. 

*****  .TEMP  CARO 

GENERAL  FORM  .TEMP  TF 1 TE2  ... 

EXAMPLE  .TEMP  -55.0  25.0  125.0 

THIS  CARD  SPECIFIES  THE  TEMPERATURES  AT  WHICH  THE  CIRCUIT  IS  TC  BF 
SIMULATEO.  TE 1 , TE2,  ...  ARE  THE  DIFFERENT  TEMPERATURES,  IN  CEGRFES  C. . A 
MAXIMUM  OF  FIVE  TEMPERATURES  ARE  ALLOWED.  TEMPERATURES  LESS  THAN  -223. C DEG  C 
ARE  IGNORED. 


AD-A055  161 


UNCLASSIFIED 


STANFORD  UNIV  CALIF  STANFORD  ELECTRONICS  LABS  F/G  9/2 

TECHNIQUES  AND  APPLICATIONS  OF  COMPUTER-AIDED  CIRCUIT  SIMULATIO— ETC (U) 
FEB  7*  R M DUTTON 
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NL 


END 


END 
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.OUTPUT  CARO 

GFNER AL  FCK*  .OUTPUT  VXXXXXX  N*  N— 

.IJUTPUT  IXXXXXX  VYYYYYY 
.OUTPUT  NOISE 

EXAMPlf  S .OUTPUT  VMIXER  13  27 

.OUTPUT  teASEl  V17 

THIS  CARO  DEFINES  AN  OUTPUT  VAR1ABUF.  FOP  VOLTAGE  OUTPUTS,  THE  NAME  MUST 
BEGIN  WITH  \ V,  4NO  N«-  ANC  N-  A RF  THE  POSITIVE  ANC  NEGATIVE  NDOF  CF  THF  OUTPUT 
VOLTAGE.  FOR  CURRENT  OUTPUTS,  THE  OUTPUT  NAME  MUST  BEGIN  WITH  AN  I,  ANC  VYYYYYY 
IS  THF  NAME  OF  THF  INDEPENDENT  VOLTAGE  SOURCE  THAT  THF  CURRENT  IS  FLOWING  IN. 
POSIT! VF  CURRENT  FLOWS  FROM  THE  POSITIVE  NOCE,  THROUGH  THE  SOURCE,  TO  THE 
NEGATIVE  NJOE.  THE  OUTPUT  VARIABLE  NAME  NOISE  IS  RESERVEO  FOR  THF  NOISF 
ANALYSIS,  AND  THE  OUTPUT  NOISE  ANO  EQUIVALENT  INPUT  NOISE  CAN  BE  PRINTED  ANO 
PLCTTEC  IN  THF  SAME  FASHION  AS  OTHER  OUTPUT  VARIA6LCS. 


OUTPUTS  CAN  BE  PRINTED  IN  TABULAR  FORM  OR  PLOTTED  AS  LINE  PRINTER  PLOTS. 
THERE  ARE  FIGHT  DIFFERENT  OPTIONS  WHICH  CAN  BE  PRINTED  ANO/OP  PLOTTEO: 


CC  OC  TRANSFER  CURVE  OUTPUT 

Tr  TRANSIENT  ANALYSIS  OUTPUT 

RF  AC  ANALYSIS  OUTPUT,  RFAL  PART 

I'M  AC  ANALYSIS  OUTPUT,  IMAGINARY  PART 

M4  AC  ANALYSIS  OUTPUT,  MAGNITUDE 

PH  AC  ANALYSIS  OUTPUT,  PHASF 

ou  noise  analysis  output,  total  output  noise  voltage 

IN  NOISE  ANALYSIS  OUTPUT,  EQUIVALENT  INPUT  NOISE 

AN  OUTPUT  GAN  BF  PRINTED  OR  PLOTTED  BY  APPENDING  THE  LETTERS  PRINT  OP  PLOT, 
FOLLOWED  BY  ANY  COMBINATION  OF  THE  EIGHT  OUTPUT  OPTIONS,  TO  THF  .OUTPUT  CARO: 


EXAMPLES  .OUTPUT  V13  13  0 PRINT  MA  DC  TR 

.OUTPUT  l I N VIN  PRINT  PH  RE  OC 
.OUTPUT  VOUT  17  2 PLOT  MA  TR  PRINT  OC 
.OUTPUT  112  V13  PLOT  PH  OC 


.OUTPUT  VTHREF  3 0 PRINT  OC  PLCT  TRAN 
.OUTPUT  NOISF  PRINT  IN  PLOT  OU 

THF  PROGPAM  WILL  AUTOMATICALLY  DETERMINE  THE  MINIMUM  AND  MAXIMLM  VALUES  OF 
THE  OUTPUT  VARIABLE  AND  SCALE  THE  PLOT  TO  FIT  THESE  LIMITS.  THF  AUTOMATIC 
SCALING  FEATURE  CAN  BE  OVER IDCEN  BY  SPECIFYING  PLCT  LIMITS  AFTER  THf  OUTPUT 
OPTION.  THE  PLOT  LIMITS  APPLY  ONLY  TO  THF  OPTION  THAT  THFY  FOLLOW. 

EXAMPLE  .OUTPUT  VI2  12  0 PLOT  MA  PH  -20  30  TR  0 5 

IN  THIS  FXAMPLF,  T EE  PROGRAM  WILL  DETERMINE  LIMITS  FOP  THE  MAGNITUDE  PLOT, 
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*****  .OC  CAPO 

gfneeal  form  .nc  op  output  input  tc  elname  vstart  vstop  vincr 

EXAMPLES  .OC  CP 

.OC  TC  VIN  0 5 0.5 

.DC  OP  VOUT  VIN  TC  VIN  0 5 0.5 

FOR  THE  SMALL  SIGNAL  TRANSFER  FUNCTION,  OUTPUT  IS  THE  OUTPUT 
VARIABLE  ANO  INPUT  IS  THE  INPUT  SOURCE.  THE  PROGRAM  MILL  CCMPUTE  THE  OC  SMALL 
SIGNAL  VALUE  OF  THE  TRANSFER  FUNCTION  ( OUTPUT/ 1 NPLT ) , INPUT  IMPECANCE,  ANO 
OUTPUT  IMPECANCE.  IF  THE  TRANSFER  FUNCTION  VALUE  IS  NOT  DESIRFO,  OMIT  THE 
OLTPLT  ANO  INPUT  SPECIFICATIONS.  IF  THE  OC  OPERATING  POINT  IS  NOT  CESIREO, 

CM  I T The  LETTERS  OP.  HOWEVER,  a OC  OPERATING  POINT  WILL  ALWAYS  BE  COMPUTEO 
PRIOR  TC  AN  AC  SMALL  SIGNAL  ANALYSIS  OR  A TRANSIENT  ANALYSIS. 

fcr  transfer  curves,  elname  is  the  name  of  the  variable  source,  vstart  is 

THE  STARTING  SOURCE  VALUE,  VSTOP  IS  THE  F INAL  SOURCE  VALUE,  ANO  VINCR  IS  THE 
INCREMENT.  THF  TOTAL  NUMBER  CF  PCINTS  TO  BE  COMPUTEC  CANNOT  EXCFED  101.  IF  A 
TRANSFER  CURVE  IS  NCT  OESIREO,  CMIT  THE  LETTERS  TC  ANO  THE  TRANSFER  CURVE 
PARAMETERS. 


*****  .AC  CARO 
Y 

GENERAL  FORM  .AC  OEC  NO  ESTART  ESTOP  NOISE  CUTPUT  INPUT  NUMS 

.AC  OCT  NO  FSTART  ESTOP  NOISE  OUTPUT  INPUT  NUMS 

.AC  LIN  NP  FSTART  ESTOP  NOISE  CUTPUT  INPUT  NUMS 

EXAMPLES  .AC  OEC  10  l 1CKHZ 

.AC  OEC  20  l 1 QOKHZ  NOISE  VOUT  VIN  10 

.AC  OEC  10  l 100MEG  NOISE  VOUT  V21- 


OEC  STANDS  FOR  OECAOE  VARIATION,  ANO  NO  IS  THE  NUMBER  OF  PCINTS  PER 
OECAOE.  OCT  STANDS  FOR  OCTAVE  VARIATION,  ANO  NO  IS  THE  NUMBER  CF  PCINTS  PER 
OCTAVE.  LIN  STANOS  FOR  LINEAR,  ANO  NP  IS  THE  NUMBER  OF  POINTS.  FSTART  IS  THE 
STARTING  FREQUENCY,  ANO  ESTOP  IS  THE  FINAL  FREQUENCY.  THE  TOTAL  NUMBER  OF 
FREQUENCY  POINTS  TO  BE  COMPUTED  CANNOT  EXCEED  101. 


FCR  NOISE  ANALYSIS,  OUTPUT  IS  THE  NAME  OF  A VOLTAGE  OUTPUT  VARIABLE.  THIS 
CUTPUT,  WHICH  MUST  BE  A VOLTAGE,  WILL  BE  USED  AS  THE  SUMMING  »CINT.  INPUT  IS 
THF  NAME  OF  AN  INDEPENDENT  VOLTAGE  OR  CURRENT  SOURCE.  THE  TOTAL  OUTPUT  NOISE  IS 
OIVIDFD  BY  THF  TRANSFER  FUNCTION  IOUTPUT/ INPUT ) TO  OBTAIN  THE  EQUIVALENT  INPUT 
NOISE  LEVEL.  NUMS  IS  THE  SUMMARY  INTERVAL.  AT  EVERY  NUMS  FREQUENCY  POINTS,  THE 
INOIVICUAL  CONTRIBUTIONS  CF  EACH  ELEMENT  ARE  PRINTED  OUT.  IF  NUMS  IS  OMITTEO 
OR  SET  TO  ZERO,  NO  SUMMARY  PRINTOUT  WILL  OCCUR.  FOR  REASONS  CF  REOUCING 
PRINTOUT,  NUMS  SHOULO  BE  AS  LARGE  AS  POSSIBLE.  IF  THF  NOISE  ANALYSIS  IS  NOT 
OFSIRFC,  OMIT  THE  LETTERS  NOISE  ANO  THE  NOISE  ANALYSIS  SPECIFICATIONS. 


.... 
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*•*«*  .TRAN  CARO 

general  FCRM 

.TRAN 

TSTEP  TSTOP  TSTART  FOUR  CUT  PUT 

EXAMPLES 

.TRAN 

INS 

1U0NS 

.TRAN 

INS 

1CC0NS  500NS 

. TRAN 

INS 

ICONS  FOUR  VOUT  100MEG 

TSTEP  IS  THE  PRINTING  INCREMENT  BETWEEN  TIMEPGINTS,  TSTOP  IS  TEE  FINAL 
TIMCPOINT,  ANO  TSTART  (S  THE  INITIAL  T I MFPO I NT.  IF  T ST  ART  IS  CMITTED,  IT  IS 
ASSURER  TC  RF  ZERO.  THE  TRANSIENT  ANALYSIS  ALWAYS  BEGINS  AT  TIME  ZERO.  IN  THF 
INTERVAL  (ZFRD,  TSTART),  Tl-6  CIRCUIT  IS  ANALYZED  (TO  REACH  A STEACY  STATE),  BUT 
NO  OUTPUTS  ARE  STORED.  IN  THE  INTERVAL  (TSTART,  TSTTP),  THE  CIRCUIT  IS 
ANALYZFC  ANO  OUTPUTS  ARE  STCRFC.  THE  NUMBER  OF  TIMFPOINTS  IN  THE  INTERVAL 
(ZERO,  ISTOP)  CANNOT  EXCEFO  1001,  ANO  THE  NUMBER  OF  TIMFPOINTS  IN  T(-g 
INTCRVAL  (TSTART,  TSTOP)  CANNOT  EXCEED  101. 

FOR  FOURIER  ANALYSIS,  OUTPUT  IS  THE  OUTPUT  VARIABLE  ANO  FRFU  IS  THE 
FUNT AMFNT AL  FREQUENCY.  TEE  FCURIER  ANALYSIS  IS  PERFORMEO  OVER  THE  INTERVAL 
( TSTOP-PERICO, TSTOP) , WHERE  TSTOP  IS  THE  FINAL  TIME  SPECIFIED,  ANO  PFRIOO  IS 
CNF  PER  ICC  OF  THE  FUNDAMENTAL  FRFCUENCY.  THF  DC  CCMPCNENT  ANO  TEE  FIRST  NINE 
COMPONENTS  ARE  DETERMINED.  FOR  MAXIMUM  ACCURACY,  THE  NUMBER  CF  PER IOOS  IN  THF 
IMTFRVAL  (TSTART, TSTOP)  SHOULC  RE  AS  SMALL  AS  POSSIBLE  (BUT  NEVER  LESS  THAN 
ONE).  THIS  INSURES  THAT  THE  NUMBER  OF  TIMEPOINTS  IN  CNE  FU  NOAMF  N TA  L IS  AS  LARGE 
AS  POSSIBLE.  IF  THE  FOURIER  ANALYSIS  IS  NOT  DESIRED,  OMIT  THE  LE  TT  E°  S FOUR  ANO 
THE  FOURIER  SPECIFICATIONS. 

FOR  SOMF  PROBLEMS,  TO  AVOID  NUMERICAL  INSTABILITY  IN  THE  INTEGRATION 
ALGORITHM,  IT  MAY  BF  NECESSARY  TO  SPECIFY  AN  INTERNAL  TIME  STEP  WHICH  IS  SMALLER 
THAN  THE  PRINTING  INCREMENT  (TSTEP).  EXAMPLES  OF  THIS  TYPE  OF  PROBLEM  ARE 
ASTABLE  MUL TI VI3RAT0RS,  SWEEP  CIRCUITS,  AND  OTHER  HIGHLY  NONLINFAR  CIRCUITS 
WHICH  HAVE  WIDCLY  SEPARATED  TIME  CONSTANTS.  SPICE  ALLOWS  THE  USER  TO  SEGMENT 
THE  TIME  INTERVAL  INTO  FRC M ONE  TO  FIVE  SUBINTERVALS  ANO  SPECIFY  A CIFFERENT 
T I MF  STEP  FOR  EACH  SUBINTERVAL.  THE  INTERNAL  TIME  STEPS  AND  SUB l N T ER VA  L 
ENCFCINTS  ARE  SPECIFIED  AFTFR  THE  STARTING  TIME  (TSTART)  ANO  BEFORF  THE  FOURIER 
ANALYST  S OPTIONS: 

GENFoau  FORM  .TRAN  TSTEP  TSTOP  TSTART  01  El  D2  E2  ...  D5  F5  FOLR  CLTPUT  FREC 

EXAMPLE  .TRAN  INS  100NS  0 O.INS  IONS  0.5NS  100NS 

01  IS  THE  FIRST  INTERNAL  T IMESTEP  ANO  El  IS  THE  ENOPOINT  OF  THE  FIRST 
SUeiNTERVAL.  02  IS  THE  SECCNO  INTERNAL  TIMESTEP  ANO  F2  IS  THE  =NOPOINT  OF  THF 
SECCND  SUBINTERVAL,  ANO  SO  CN.  IN  THIS  EXAMPLE,  THE  PROGRAM  wILL  USF  AN 
INTERNAL  T!Mg  STEP  OF  O.INS  FC«  THE  INTERVAL  (0.10NS)  ANO  AN  INTERNAL  TIME  STEP 
OF  C.5NS  FOR  THF  INTFRVAL  ( 10  N S , 1 OONS  ) . OUTPUT  IS  STUL  STCRFD  FVERY  INS.  THE 
TOTAL  NUMBER  OF  TIMEPOINTS  TO  BE  COMPUTED  CANNOT  EXCEED  1001. 

F X A*»PL  F .TRAN  ILS  100US  0 O.IUS  100US 

in  this  example,  the  prcgram  will  use  an  internal  time  step  of  c.ils  over  the 

ENTIRE  TRANSIENT  INTERVAL  BUT  WILL  STORE  OUTPUT  ONLY  AT  1US  INTERVALS.  HENCE, 
THE  PROGRAM  STORES  ANO  OUTPUTS  EVERY  TFNTH  TIMFPOINT. 
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EXAMPLF  DATA  DECKS 


THE  FOLLOWING  CECK  CETERMINFS  THE  CC  OPERATING  PCINT  ANP  SPALL  SIGNAL 
TRANSFER  FUNCTION  OF  4 SIMPLF  DIFFERENTIAL  PAIR. 

SIMPLE  DIFFFRPNTIAL  PAIH 

VCC  7 C PC  12 

VF?  8 0 PC  -12 

VIN  1 0 

RSI  1 2 IK 

RS2  fc  0 IK 

01  1 2 A modi 

02  5 6 A MOCl 
RCl  7 3 10K 
RC2  7 5 10K 
RE  A 8 10K 

.MODEL  MODI  NPN  9F=50  V A=50  IS*l.0E-l2  RB*lOO 
.OLT  VCLT  5 0 
.CC  CP  VOLT  VIN 
.FNO 


THF  FOLLOWING  DECK  DETERMINES  THE  DC  TRANSFER  CURVE  ANC  THE  TRANSIENT 
PULSE  RESPONSF  OF  4 SIMPLE  RTL  INV6RTFR.  THE  INPUT  IS  A PULSE  FRCM  0 TO  5 
WITH  OELAY « RISE,  AND  FALL  TI*»ES  CF  2NS  ANO  A PULSF  WIDTH  OF  30NS.  THE 
TRANSIENT  INTFRVAL  IS  0 TO  100NS  IN  INS  STEPS. 

SIMPLE  RTL  INVFRTFR 
VCC  A 0 OC  5 

VIN  1 0 PULSF  0 5 2NS  2 N S 2NS  30NS 
R8  l 2 10K 
01  3 2 0 01 
RC  A 5 IK 

.OLTPUT  VC  3 0 PRINT  DC  PLOT  TR  0 5 

.MODEL  01  NPN  3 F *20  R8* 1 00 . TF*0 . INS  CJC*2PF 

.CC  TC  VIN  C 5 C.l 

.TRAN  INS  100NS 

.END 

THE  FOLLOWING  CECK  DETERMINES  THE  AC  SMALL  SIGNAL  RESPONSE  OF  A ONE 
TRANSISTOR  AMPLIFIER  OVER  THE  FREQUENCY  RANGE  OF  IHZ  TO  100MEGHZ. 

CNF  transistor  amplifier 

VCC  5 0 DC  12 

VEE  6 0 DC  -12 

VIN  l 0 AC  l 

"S  1 2 IK 

01  2 2 A X 3 3 

RC  5 3 500 

RE  A 6 IK 

CBYPASS  A 0 IUFO 

.OUT  V3  3 0 PLOT  MA  PH 

•AC  DEC  10  IHZ  100MEGHZ 

.MCCEL  X33  NPN  BF«30  RB*50  VA*20 
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External  Models  In  SPICE 


SPICE  allows  external,  user-defined  models.  This  feature  is  particularly 
convenient  when  a large  circuit  includes  several  identical  subcircuics  such  as 
operational  amplifiers  or  logic  gates.  The  external  model  card  is  defined  by  a 
..MODEL  card,  a set  of  element  cards  (which  may  include  any  legal  SPICE  element 
except  reference  to  another  external  model),  and  a .FINIS  card.  The  .MODEL 
card  format  is 


.MODEL  model  name  X node  1 node  2 ...  node  N 

The  .MODEL  card  Introduces  subsequent  cards  as  a definition  of  an  external  model. 

All  subsequent  cards  up  to  the  next  .FINIS  card  are  treated  as  a single  definition. 
These  cards  are  conventional  element  and  device  cards  and  refer  to  the  node  numbers 
that  appear  in  the  introducing  .MODEL  card.  Note  that  with  the  exception  of  node 
zero,  the  node  numbers  that  are  used  in  the  definition  of  an  external  model  are 
dummy  node  numbers  and  may  be  given  the  same  number  as  used  within  the  circuit  to 
be  analysed. 

Reference  to  an  external  model  is  similar  to  reference  to  a built-in  model. 
However,  the  names  of  devices  described  by  an  external  model  must  begin  with  the 
letter  X (just  as  transistors  device  names  must  start  with  the  letter  Q) . Thus 
one  writes 

X device  name  nl  n2  . . . nn  model  name 

The  nodes  nl  thru  nn  are  actual  circuit  nodes.  They  are  paired  with  the  node  numbers 
of  the  external  .MODEL  definition  on  the  basis  of  position  alone.  The  first  node 
in  the  model  card  is  paired  with  the  first  node  in  the  device  card,  the  second  with 
the  second  and  so  on. 


»«c  » a oc  l.4» 

»l'l  1 0 Kilt  w 1.4*  JM  CVS  CM  VMS 

• T44*  ms  1SCSS 

.jut  **«o  nvr  r**i 

• 

«N»I  12  5 IMKMfH 
2)5  invcOff. 
u.v)  )55  hvnn 
• 

• A 10  AC 

n jo  20  a »ec 

*4  U 20  *50 
<C  a *3  JC  (AC 

• *ica  MC  **1  *M)4  MM  «0MJ*  «CM0O 

• M«t*$  MmjnS  CJCMM  CJ'«2I»* 

••TNI S 


f 

t 

\ 


